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We develop an inertial coupling method for modeling the dynamics of point-like "blob" 
particles immersed in an incompressible fluid, generalizing previous work for compressible 
fluids [F. Balboa Usabiaga, I. Pagonabarraga, and R. Delgado-Buscalioni, J. Comp. Phys., 
235:701-722, 2013]. The coupling consistently includes excess (positive or negative) inertia 
of the particles relative to the displaced fluid, and accounts for thermal fluctuations in the 
fluid momentum equation. The coupling between the fluid and the blob is based on a no- 
slip constraint equating the particle velocity with the local average of the fluid velocity and 
conserves momentum and energy. We demonstrate that the formulation obeys a fluctuation- 
dissipation balance, owing to the non-dissipative nature of the no-slip coupling. We develop 
a spatio-temporal discretization that preserves, as best as possible, these properties of the 
continuum formulation. In the spatial discretization, the local averaging and spreading 
operations are accomplished using compact kernels commonly used in immersed boundary 
methods. We find that the special properties of these kernels make the discrete blob a particle 
with surprisingly physically-consistent volume, mass, and hydrodynamic properties. We 
develop a second-order semi-implicit temporal integrator that maintains discrete fluctuation- 
dissipation balance, and is not limited in stability by viscosity. Furthermore, the temporal 
scheme requires only constant-coefficient Poisson and Helmholtz linear solvers, enabling a 
very efficient and simple FFT-based implementation on GPUs. We numerically investigate 
the performance of the method on several standard test problems. In the deterministic 
setting, we find the blob to be a remarkably robust approximation to a rigid sphere, at both 
low and high Reynolds numbers. In the stochastic setting, we study in detail the short and 
long-time behavior of the velocity autocorrelation function and the mean square displacement 
of a freely diffusing particle. We find a surprising deviation from the standard Stokcs- 
Einstein prediction at small Schmidt numbers, owing to a non-averaging of the fluctuating 
contributions to the particle mobility. The proposed inertial coupling method provides a 
low-cost coarse-grained (minimal resolution) model of particulate flows over a wide range of 
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time-scales ranging from Brownian to convection-driven motion. 

I. INTRODUCTION 

The dynamics of small particles immersed in a fluid is key to many applications involving dis- 
parate length and time scales pQ: from the dynamics of millimeter particles (dust) in turbulent 
flow, to multiphase flow with micron and nanoscopic colloidal molecules in quiescent, laminar [Ji- 
ll] , or turbulent regimes [H [6] . In many engineering applications colloidal particles are exposed to 
disparate dynamic regimes coexisting in different subdomains of the same chamber [7J. Such pro- 
cesses demand fast computational methods able to efficiently resolve the motion of many (O (10 5 )) 
colloidal particles in quite different dynamics ranging from diffusive to inertial dynamics. Such 
scenarios are paradigmatic of what one might call multi-regime systems. 

A group of methods such as smooth particle hydrodynamics (SPH) [8j, smoothed dissipative 
particle dynamics [9], and stochastic rotation dynamics (SRD) [2] resolve both the particle and fluid 
phase using similar discrete Lagrangian descriptions, and as such, seem to be natural candidates 
to become multi-regime solvers [8j [10] . Particle-particle methods allow for an easy treatment of 
complex boundary conditions, and offer a natural way to couple moving boundaries or immersed 
particles to the fluid. However, particle-particle methods have important drawbacks when com- 
pared with standard solvers for discretized Computational Fluid Dynamics (CFD). In particular, 
they offer limited control over the fluid properties and require relatively small time steps com- 
pared to, for example, semi-implicit CFD schemes. Moreover, they cannot be adapted to efficiently 
treat the natural time scales governing the different dynamical regimes (e.g., the Brownian or 
overdamped limit). Similar advantages and drawbacks also apply to the lattice Boltzmann (LB) 
method [3J, although the LB approach has proven to be a rather flexible framework [llj. 

Many other approaches use CFD for the solvent flow and couple its dynamics with that of 
the immersed particles. In the realm of CFD one can still distinguish two large subgroups of 
methods. The first group of methods involves a Lagrangian description of the computational mesh 
which self-adapts to follow the particle |12j . The second group uses a fixed (Eulerian) grid and 
requires converting the particle boundary conditions into body forces or some interaction equations 
[TTl HU E]. The present work focuses on this second group, sometimes called mixed Eulerian- 
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Lagrangian methods. These schemes are particularly suited to attack the "multi-regime" problem, 
because they are faster, more flexible and can work with minimal resolution models (pointwise 
particle descriptions) . 

In their seminal work, Maxey and Riley [15] decomposed the fluid velocity as v(r, f) = v (r,i) + 
Vi(r, t), where v is the undisturbed flow (which would result if the boundary conditions at the 
particle surface were not applied), and v x is the perturbative component created by the fluid- 
particle interaction. In the bulk flow, convection (advection) becomes relevant for Re P = v Lp/rj > 1; 
where the fluid Reynolds number Re P is defined in terms of the typical flow speed v , the fluid 
density p, the dynamic viscosity rj = pis, and a characteristic length L for velocity variation in the 
flow. The fluid force on the particle arises from the local fluid inertia (proportional to the local 
material derivative of v ) and also from the local stress created by the particle disturbance. The 
relaxational part of the particle inertia is a consequence of its mass resistance to instantaneously 
follow the velocity of the surrounding fluid; the fluid drag damps the particle velocity towards 
the local fluid velocity within an inertial time r P ~ (p P — p)R 2 /r) which increases with the density 
contrast p P — p and with the particle radius R. 

By contrast, convective inertia arises from non- linear interactions between the particle dynamics 
and its perturbative flow [16]. The particle Reynolds number Re P = 2wR/v, defined with the 
particle- fluid relative speed w [E] , determines the relative strength of advection by the perturbative 
flow relative to viscous dissipation. The importance of convective inertia is indicated by the ratio 
Re P (i?/L) 2 between the characteristic times associated with Stokes drag and convection 1 1 5 1. fTGlj . At 
finite values of the non-dimensional groups Re P and Re F (R/L) 2 inertia effects due to particle mass 
and particle size are not interchangeable anymore, especially in the turbulent regime [SJ[T7]. Non- 
linear interaction between particle advection and thermal fluctuations are also possible at small 
Reynolds number. Some examples are the change in the mobility of colloidal particles R ~ 10~ [5 ~ s] m 
over the Stokes limit at low values of the Schmidt number (typical of aerosol) |18[ 119] . and inertial 
effects in directional locking, a process to to separate nanoparticles at very small Re P [20 . 

Computational approaches are usually tailored to tackle some specific dynamical regime and 
they can be naturally classified according to the range of Re F , Re P and R/L they can be safely 
applied to. In the creeping flow limit, Re P — > and Re P — > 0, the perturbative flow vi has a 
negligible effect on the unperturbed field, which is a priori fixed. The perturbative field created by 
a collection of particles is the linear superposition of the Stokes fields and it determines the multi- 
body hydrodynamic forces on the particle ensemble. Analytical expressions for these forces are 
embedded in the mobility matrix of Brownian hydrodynamics (BD) [21|. 122] and Stokesian dynamics 
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(SD) [23J. In addition to the stokeslet (monopole) terms, in SD one can include higher terms of 
the multipole expansion of the perturbative stress pQ . The zero- Reynolds regime resolves the long- 
time diffusive (Smoluchowski) limit of colloidal motion, in which fluctuations make an important 
(O(l)) contribution. Direct implementation of the fluctuation dissipation (FD) relation between 
the friction and noise matrices requires 0(N 3 ) operations, where N is the number of particles. 
Sophisticated and technically-complex techniques such as the accelerated Stokesian dynamics [24] , 
and the general geometry Ewald-like method for confined geometries [25] , reduce the large raw cost 
to 0(N In TV) operations, albeit with large multiplicative prefactors. 

As an alternative to BD and SD methods, two-way coupling algorithms using a Stokes frictional 
force were developed for mixed Eulerian-Lagrangian dynamics [3] [26] [27] . The idea is to deploy 
a relative simple and efficient fluid solver to explicitly resolve the perturbative flow responsible 
for the hydrodynamic coupling between particles. The particle cost is dominated by neighbour 
searching and scales (almost) linearly with N, while the (added) fluid solver cost scales like the 
system volume. The Eulerian-Lagrangian mixed approach permits to work at finite Re P . However, 
the Stokes (i.e. frictional) coupling assumption limits the scheme to Re P < 1 and only resolves 
far-field hydrodynamics (R/L < l). The Stokes coupling consistently neglects convective inertia 
and only includes some relaxational particle inertia, and the finite particle response time t> is 
introduced via a phenomenological friction coefficient. Frictional coupling is obviously dissipative 
and requires introducing an additional noise term in the particle equation, different from that of 
the fluctuating fluid stress tensor [3] [27] . 

Other methods for finite Re P have been restricted to Re P = 0, where the particle inertia is absent. 
Two relevant examples are the stochastic Immersed Boundary method (IBM) [271 [2"B] commonly 
used for fluid-structure interaction at R/L = 0(1), and the Force Coupling method (FCM) [15] 125], 
where each particle is represented by a low-order expansion of force multipoles (R/L < 1). An 
extension of the FCM that consistently includes thermal fluctuations is, to our knowledge, still 
lacking. For Re P = the relative fluid-particle acceleration is zero and the particle velocity just 
follows the local fluid velocity. The hydrodynamic force due to the particle-fluid interaction is then 
equal to the total force exerted on the particle by sources other than the fluid. This permits a 
fluid-only formulation whereby the net non-hydrodynamic particle force is spread from the particle 
to the surrounding fluid using some compact kernel. This important spreading operation differs 
substantially from method to method. In FCM two different Gaussian kernels are used to spread 
the force monopole and force dipole moments (stresslet); their widths are fitted in the continuum 
model to recover the Stokes drag and linear Faxen terms [29] . By contrast, the IBM kernels are 
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specifically designed to minimize the effects of the discrete Eulerian mesh in the spreading of 
Lagrangian point forces (monopole terms) |30| . 

For Re P <C 1 and R/L <c 1 advection of the perturbative flow can be neglected leading to the 
(analytically solvable) unsteady Stokes equation for the perturbative field [15]. The fluid-particle 
force can be expressed as some function of the relative velocity field u— v interpolated at the particle 
site. This forms the basis of one-way- coupling schemes for point-particle dynamics frequently used 
in turbulence research Re P ^> 1 |16j . Generalizations to Re P ~ 1 have been also derived (see e.g. 
|31j ) but, even in the simpler Re P <c 1 limit, the evaluation of the fluid-particle force involves 
cumbersome expressions which require interpolations of the displaced fluid acceleration (the added 
mass effect), second order spatial derivatives of v (Faxen terms), and time-convolved integrals 
which recast the history of vorticity diffusion around the particle (Basset memory). For a sphere 
moving with velocity u at Re P <C 1, the leading term is the steady Stokes force F S tokes = QitnR(xi— v ), 
which, due to its simple form, has been overused in two-way point-particle approximations of 
turbulent (Re F ^> l) and pre-turbulent regimes [3[T7]. Although the point-particle approach can 
probably describe the relaxational inertia of very small {R/L <g; 1) heavy particles in a light fluid 
(e.g. aerosol), it has the serious limitation of neglecting the convective inertia arising from the 
particle finite size |32j . Even at low Re P , convection of perturbative flow is known to alter the 
Basset memory and the long time particle dynamics [16J. And vice versa, recent works show that 
micron-size particles can alter the turbulent spectra at moderate Re P and R ~ L (with L the 
Kolmogorov length) [5j [6] due to energy dissipation and vorticity production in the particles wake 
[32]. 

Several Eulerian-Lagrangian methods have appeared in recent years to allow for a fully consistent 
treatment of the coupled particle and fluid inertia. A key issue is the spatial resolution of the 
particle. In the "direct forcing" method [33], and related extensions to fluctuating hydrodynamics 
[3"H 135] , the fluid force on the particle is obtained by imposing the no-slip constraint on a 
well-resolved particle surface (and perhaps also the interior of the particle) . High spatial resolution 
requires substantial computational effort; the largest simulations so far reached O(10 s ) particles 
|33tl36|. The smoothed particle method (SPM) |13tll4| works with a mixed (particle-fluid) velocity 
field constructed with a smooth characteristic function which discriminates particle and fluid cells. 
This permits an intermediate resolution with a typical particle radius R ~ 5h (here h is the mesh 
size) requiring O(10 3 ) fluid cells per particle. 

These fully or partially resolved methods are still far from a point-particle approach which 
can require as few as 13 cells to perform a fourth order orthogonal Lagrangian interpolation [37J. 
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Pointwise or "blob" particle descriptions (discussed in more detail in the next section) offer a way to 
explore finite particle effects at moderate computational cost. In a preceding paper |38| . some of us 
proposed an inertial coupling method that directly couples a compressible finite- volume fluctuating 
hydrodynamic solver [39j to blob particles. A distinguishing feature of the coupling methodology 
is that it includes the effect of the particle and fluid inertia in the dynamics, while still consistently 
including thermal fluctuations even in non-trivial geometries. It was numerically demonstrated 
that the inertial coupling method can reproduce ultrasound forces on colloidal particles, taking 
place at much faster rates than viscous friction |40j . 

In the previous work [38J, a compressible solver was used because one of the focus applications 
was the interaction between ultrasound and colloidal particles |40| . In many applications sonic 
effects can be ignored and the essential hydrodynamic interactions can be captured by using the 
isothermal incompressible Navier-Stokes equations instead of the compressible equations. This 
eliminates the fast sound waves and allows for a much larger time step size in the fluid solver. Here 
we develop an Inertial Coupling Method (ICM) that directly couples an incompressible finite- volume 
solver for the fluctuating Navier-Stokes equations (39] with suspended particles, which do not 
necessarily have the same density as the fluid. We demonstrate that the coupling obeys a continuum 
and a discrete fluctuation-dissipation balance and study the performance of our algorithm. The 
ICM is a coarse-grained model for particle hydrodynamics which aims to capture hydrodynamic 
effects (unsteady forcing, viscous friction and advection) over a broad range of time scales and Re P : 
from Brownian motion to convection-driven regimes. In the ICM, the coupling between the particle 
and the fluid is not assumed to have any functional form (e.g. Stokes drag) but naturally arises 
from the no-slip constraint averaged over the particle (or "blob") domain. The present results (see 
also Ref. [38]) indicate that this type of (non- linear) coupling permits to take into account both 
fluid and particle inertia beyond the Stokes limit, where advective interactions take place. 

In the remainder of this Introduction we introduce some notation and fundamental concepts. 
In Section [IT] we discuss the continuum equations of the incompressible inertial coupling method. 
We present both a constrained and a constraint-free formulation, demonstrate momentum and 
energy conservation, as well as fluctuation-dissipation balance, and discuss the Brownian (diffusive) 



limit. In Section III we present a second-order semi- implicit spatio-temporal discretization of the 



continuum equations, and summarize our algorithm. In Section IV we test and apply the algorithm 
to a collection of standard test problems. We begin by demonstrating second-order temporal 
accuracy in the deterministic setting. We also study in detail the short and long-time behavior of 
the velocity autocorrelation function and the mean square displacement of a freely diffusing particle. 
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We then demonstrate that the algorithm correctly reproduces static equilibrium properties such as 
the radial distribution function in a suspension of soft spheres, and also reproduces known features 
of pairwise hydrodynamic interactions between two particles. Finally, we study the behavior of 
the blob particle in high Reynolds number flow. In Section [V] we offer some conclusions and 
thoughts on possible extensions of the method and algorithm. Several more technical calculations 
and algorithmic details are presented in a collection of Appendices. 

A. Notation and Basic Concepts 

In the beginning, we focus on the continuum formulation of the fluid-particle coupling. However, 
it is important to point out that most of the notation and conclusions can directly be adopted in 
the discrete formulation by simply replacing spatial integrals with sums over grid points. We will 
return to the spatially-discrete formulation in Section [TTTj 

Let us consider a particle of physical mass m and size (e.g., radius) a immersed in a fluid with 
density p. In real problems there will be many particles i = 1, . . . , N p that interact with each other, 
for example, in microfluidic applications involving polymers each particle could represent a bead 
in a bead-spring or bead-link polymer model [22J. Unless otherwise indicated it is straightforward 
to extend the proposed formulation to a collection of interacting particles by simply adding a 
summation over the different particles. Therefore, for simplicity of notation, we will typically focus 
on a single particle and omit the particle index. 

The position of the particle is denoted with q(t) and its velocity with u = q. The shape of the 
particle and its effective interaction with the fluid is captured through a smooth kernel function 
5 a (Ar) that integrates to unity and whose support is localized in a region of size a. For example, one 
may choose any "bell-shaped" function with half-width of order a. In immersed-boundary methods 
|30| . the kernel function S a is considered to be an approximation of the Dirac delta function of 
purely numerical origin. By contrast, in the force-coupling method [29|, the shape of the kernel 
function is related to the physical size and properties of the actual particle. We adopt an approach 
that is intermediate between these two extremes and choose the shape of the function based on 
numerical considerations, but relate its shape to the physical properties of the particle. 

The interaction between the fluid and particle is mediated via the kernel function through two 
crucial local operations. The local averaging linear operator J(q) averages the fluid velocity inside 
the particle to estimate a local fluid velocity 



s 



v q (t) = Jv(r,t) = I S a (q - r)v(r,t) dr 



The reverse of local averaging is accomplished using the local spreading linear operator S(q) which 
takes a force F applied to the particle and spreads it over the extent of the kernel function to 
return a smooth force density field, 



For notational simplicity we will slightly abuse notation and assume that the local spreading and 
interpolation operators can be applied to a scalar, a vector, or a tensor field, with the interpretation 
that the same local averaging or spreading operation is applied to each component independently. 
This sort of block-diagonal form of the spreading and interpolation operators is not strictly required 
for the mathematical formulation |27j , but applies to the specific Peskin forms of the operators we 
use in practice [30] . 

The physical volume of the particle AV is related to the shape and width of the kernel function 

via 



Therefore, even though the particle is represented only by the position of its centroid, it is not 
appropriate to consider it a "point" particle. Rather, it can be thought of as a diffuse sphere that 
has some physical extent and interacts with the fluid in its interior. For lack of better terminology, 
we will refer to such a diffuse particle as a "blob". In fluctuating hydrodynamics the fluid velocity 
is a distribution and cannot be evaluated pointwise, therefore, to obtain well-defined fluctuating 
equations spatial averaging must be used and a physical volume associated to each blob. 

The fluid velocity field is denoted with v(r,t) and is assumed to extend over the whole domain, 
including the particle interior. Because fluid permeates the interior of the particle, the effective 
inertia of the particle is enlarged by pAV giving the physical mass 



where m e is the excess mass of the particle over the mass of the entrained fluid m f = pAV . In 
particular, m e = corresponds to a neutrally-buoyant particle, meaning that the inertia of the fluid 
is unchanged by the presence of the particle. It is a crucial property that AV is a constant that only 
depends on the shape of the kernel function and not on the position of the particle, and therefore 
the mass of the particle m is constant and can be given a well-defined physical interpretation. 



f(r,t) = SF(t) = F(t) S a (q-r). 




(1) 



777 = 777. e + pAV 



m e + To/, 
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Note that the local spreading operator S has dimensions of inverse volume. One could alterna- 
tively use the dimensionless operator S = SAV, with the property that JS = 1. We prefer to use 
the dimensional version because the averaging and spreading operators are adjoint, S — J* , i.e., 
the natural dot products in the particle (Lagrangian) and fluid (Eulerian) domains are related via 

m 

(Jv) ■ u — J v ■ (Su) dr — J S a (q — r) (v ■ u) dr (2) 

for any u and v. This adjoint property is crucial in maintaining energy conservation and fluctuation- 
dissipation balance and will also be preserved by the discrete local averaging and spreading oper- 
ators. 



B. Fluctuating Incompressible Navier-Stokes Equation 

In this work we assume that the fluid can be described via the fluctuating Navier-Stokes equation 
|41| . Specifically, we model the dynamics of the fluid velocity field v(r,t) assuming an isothermal 
incompressible Newtonian fluid, V • v = 0, 

p (d t v + v-Vv) = - Vtt - V • a + f = - Vtt + r,V 2 v + V • UkgTv) 112 (W + W T ) ] + f, (3) 
where the stress tensor er includes viscous and fluctuating contributions, n is the non- 
thermodynamic pressure, p is the (constant) density, rj — pv is the (constant) shear viscosity, v 
is the kinematic viscosity, and / (r,t) is an additional force density such as gravity or the force ex- 
erted by the particles on the fluid. Note that we prefer to use the standard physics notation instead 
of the differential notation more common in the mathematics literature since there is no difference 
between the Ito and Stratonovich interpretations of stochastic integrals for additive noise. 

In the momentum conservation law ([3]), the stochastic momentum flux is modeled using a 
white-noise random Gaussian tensor field W(r,t), that is, a tensor field whose components are 
independent (space-time) white noise processes, 

(W«(r 1 *)W w (r',0) = & k 6 3l )5(t-t')6(r-r'). 
The form of the stochastic forcing term ensures fluctuation-dissipation balance between the random 
forcing and the viscous dissipation and gives the correct spectrum for the thermally-induced velocity 
fluctuations. The symmetrized form of the fluctuating stress (k B Tri) 1/2 (W + W T ) mimics the 
symmetry of the viscous stress is w (Vv + V T v) ; however, for incompressible flow with constant 
viscosity one may also use the stochastic flux (2k B Tr]) 1/2 W The discretization and numerical 
solution of ([3]) is discussed in more detail in Refs. [39| EE2]. 



10 



C. No-Slip Condition 

Coupling of a continuum (fluctuating) fluid with point-like (blob) particles has been considered 
by other researchers. In particular, in Lattice-Boltzmann methods 0, H3"l B3] a Stokes frictional 
force between the particle and the fluid is postulated. Specifically, the motion of the particle is 
described by a Langevin equation in which a phenomenological Stokes frictional force between the 
particle and the fluid is postulated, proportional to the difference u — Jv between the particle and 
the locally-averaged fluid velocity. A corresponding force is added to the fluid equations to ensure 
momentum and energy conservation and fluctuation-dissipation balance in the fluid-particle system 

13 E7J. 

An important downside of the inertial Stokes coupling is the imposition of an artificial friction 
parameter and an associated delay with the response of the particle to changes in the flow. Such 
a delay is often not physically acceptable unless a very large friction constant is imposed, leading 
to numerical stiffness. Instead, following Ref. [38], we impose an instantaneous coupling between 
the fluid and the particle in the form of a no-slip constraint, 

u = q = Jv, (4) 
The no-slip condition simply states that the velocity of the particle is equal to a local average 
of the fluid velocity. This is a constraint that formally eliminates the particle velocity from the 
formulation and leaves only the fluid degrees of freedom. We now demonstrate that the imposition 
of Q leads to a physically-consistent and sensible coarse-grained model of the coupled fluid-particle 
system. Notably, our coupling conserves momentum, energy, and obeys a fluctuation-dissipation 
principle. 

The particle acceleration is 

u=j t [J{q)v] = J{d t v)+(u.^j} V , (5) 
where for our choice of interpolation operator we have the explicit form: 



u ■ -^-J ] v = 
oq 



u ■ §- q 5 a (q r) 



v (r, t) dr. 

Observe that in the limit of a "point particle", a — > 0, the kernel function approaches a Dirac delta 
function and one can identify ^ with the advective derivative, 

~ {Jv) « ~v (q(t),t) = D t v = d t v + (v ■ V) v, 
at at 

which is expected since in this limit the particle becomes a Lagrangian marker. In Ref. [45]. 
the term ^ (Jv) is replaced with the interpolated Navier-Stokes advective derivative J (D t v), thus 
avoiding the need to differentiate the kernel function. For a blob particle, however, in general, the 
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relative fluid-particle acceleration is non-zero, 

aj = Jt (J?;) ~ J ^ DtV) = ( u ' Jg J ) V JV ' Vv * ° ^ 

II. INCOMPRESSIBLE INERTIAL COUPLING METHOD 

Following the discussion in the Introduction and the derivation in Section 2 of Ref. [38] we take 
the equations of motion for a single particle coupled to a fluctuating fluid to be 

p (d t v + v ■ V«) = pD t v = -Vtt - V • a - S (q) A (7) 
m e u = F(q) + X (8) 
s.t. u = J (g) (9) 
where the fluid-particle force A is a Lagrange multiplier that enforces the constraint ([9]) and F (g) 
is the external force applied to the particle. Observe that the total particle-fluid momentum 



P 



J pv (r, t) dr 



is conserved because Newton's third law is enforced; the opposite total force is exerted on the fluid 
by the particle as is exerted on the particle by the fluid. When there are more than one particle 
one simply adds the forces from all the particles in the fluid equation. 

Note that similar equations apply for both compressible and incompressible fluids. In the 
compressible case [HH], a density equation is added to the system ( 7|8|9 ) and the pressure ir (p) 



obtained from the equation of state. In the incompressible case the divergence-free condition 
V • v = is used instead to determine the (non-thermodynamic) pressure as a Lagrange multiplier. 

For now, we will silently ignore the fact that in the fluctuating equations includes a non-smooth 
white noise component that must be handled with care, and return to a discussion of the stochastic 
equations later on. For a neutrally-buoyant particle, m e = 0, A = — F, and the fluid equation is the 
standard Navier-Stokes with the force on the particle spread back to the fluid as a force density SF 
\29\ [30] . In this case our formulation is equivalent to the Stochastic Immersed Boundary Method 
[271 [28]. 

A. Effective fluid and particle equations 



In this section we study the properties of ( 7|8|9 ) in order to better understand the physics of 



the fluid-particle coupling. Using ^ to eliminate A = m e it — F and Q to eliminate it, the fluid 
equation ([7]) becomes, 

pD t v = p{d t v + v ■ Vv) = -m e SJ (D t v) - Vtt -V ■ a -m e Saj + SF. (10) 
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This gives the effective fluid equation 



(p + m e SJ) d t v = — 



/'If V! + in, S [u - J^J 



v- Vtt- V-er + SF, (11) 



in which the effective fluid inertia is given by the operator p + m e SJ, and the kinetic stress term 
pv ■ Vi> includes an additional term due to the excess inertia of the particle. When there are many 
interacting particles one simply adds a summation over all particles in front of all terms involving 



particle quantities in (11). Note that for a neutrally-buoyant particle m e = and one obtains the 
constant-density Navier-Stokes equation with external forcing SF. 

Similarly, by eliminating A from Q we obtain the effective particle equation (see also Section 
2 of Ref. [38]), 

mil = ~AV J (Vtt + V ■ er) + F + m f a,j, (12) 

where m f = pAV is the mass of the fluid dragged with the particle. This equation makes it clear 
why m — m e + m,f has the physical interpretation of particle mass (inertia) . If the particle were a 
rigid sphere, the force exerted by the fluid on the particle would be the surface average of the stress 
tensor. It is sensible that for a blob particle this is replaced by the locally averaged divergence of 
the stress tensor (first term on right hand side). The last term in the particle equation m f aj has a 
less-clear physical interpretation and comes because the fluid is allowed to have a local acceleration 
different from the particle. It is expected that at small Reynolds numbers the velocity field will be 
smooth at the scale of the particle size and thus aj ss [36]. Nevertheless, we will retain the terms 
involving a,j to ensure a consistent formulation, see Appendix [B] 



B. Momentum Conservation 

Let us define a momentum field as the sum of the fluid momentum and the spreading of the 
particle momentum, 

p (r, t) = pv + m e Su = (p + m e SJ) v. (13) 

The total momentum is P (t) = J p(r, t) dr and therefore a local conservation law for p(r,t) implies 
conservation of the total momentum. 

By adding the fluid and particle equations 
momentum field, 

d t p = p (d t v) + m e Su + m e ( u ■ S } u 

V dq ) 

= -Vn - V • a - V • [pvv T + m e S (uu T )] + SF, (14) 



7p) together we can obtain the dynamics of the 
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where we used the fact that S depends on the difference (q — r) only, and not on q and r separately. 
In the absence of applied external forces we can write the right hand side as a divergence of a total 
stress tensor nl + rr + er kin , where the kinetic stress tensor includes a contribution from the inertia 
of the particle, 

Ckin = pvv T + m e S (im T ) . (15) 

This means that the momentum field obeys a local conservation law, as expected for short-ranged 
interactions between the particle and the fluid molecules. 



The formulation (14) is not only informative from a physical perspective, but was also found 
very useful in performing adiabatic elimination in the case of frictional coupling in Ref. [27]. We 
have also found the momentum formulation well-suited for constructing numerical discretizations 
that exactly discretely conserve momentum. 

C. Energy Conservation 



In the absence of viscous dissipation, the equations of motion ( 7|8f 9) conserve a coarse-grained 



Hamiltonian \47\ |4"8] given by the sum of potential energy and the kinetic energy of the particle 
and the fluid, 

— dr + m a — + U{q) 1 (16) 

where U (q) is the interaction potential of the particle with external sources and other particles, 
with an associated conservative force 

, . dU dH 
F ^ = -dq=-W 

For compressible flow one needs to include the (density-dependent) internal energy of the fluid in 
the Hamiltonian as well 



To demonstrate energy conservation, we calculate the rate of change 

= — F ■ u + m,,u ■ u + / pv ■ (o t v) dr 

at J 

in the absence of viscous and stochastic fluxes. Using the equations of motion ( 7[8 ) we get 

dH 

lit 



-F u + u- (F + A) - J v- (SX) dr 
— J v ■ V-7T dr — p j v ■ (v ■ Vw) dr 
(u - Jv) • A + j 7r (V • v) dr 
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where integration by parts and the adjoint property ([2]) were used for the first two terms, and a 
vector identity was used to express w-V« in terms of the vorticity V x v. The first term vanishes due 
to the no-slip constraint u = Jv. The second and third terms vanish for incompressible flow V • v, 
and the last term vanishes because of the basic properties of the cross product. This demonstrates 
that dH/dt = in the absence of viscous dissipation, that is, the non-dissipative terms in the 
equation strictly conserve the coarse-grained free energy. 



D. Pressure- Free Formulation 



The equations we wrote so far contain the V-7T term and can easily be generalized to the case of 
a compressible fluid [38]. For analysis purposes, in the incompressible case it is useful to eliminate 
the pressure from the equations using a projection operator formalism. This well-known procedure 
can be understood as follows. The fluid equation ([7]) is of the form 

d t v + p^V-K = g, V-v = 0. 
By taking the divergence of the evolution equation, we get 

dt (V • v) + p- l V 2 ir = p^VV = V • g, 
which is a Poisson equation for the pressure whose solution can be formally written as n = 
pV~ 2 (V ■ g). This means that 

d t v = g- p-'Vn = g - V [V~ 2 (V ■ g)] = Vg. 
where V is a projection operator that projects the right-hand side g or a given velocity field onto 
the space of divergence-free vector fields. Note that the boundary conditions are implicit in the 
definitions of the gradient, divergence, and Laplacian operators. For periodic boundary conditions, 
the projection can most easily be implemented using a spatial Fourier transform. Specifically, in 
Fourier space the projection operator is simply a multiplication by the dxd matrix V = I—k~ 2 (kk T ), 
where d is the dimensionality, d = 2 or d = 3, and k is the wavenumber. 

By using the projection operator, we can eliminate the pressure from the equations of motion 
(7p), to obtain 



pd f v = V [-pv ■ Vv - V • a - m e Sii + SF] . 

If we now use §5§ to eliminate u we obtain the fluid equation 

pd t v + m e VSJV (d t v) = V -pv ■ Vv - m e S ■ v - V ■ a + SF 

where we used the fact that Vv = v since V • v = 0, and we added a V in front of the second term 
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for symmetry purposes. This shows that the pressure- free form of the fluid-only equation (11) is 

p eS d t v ="Pj- P(V- V) + m e SJ\v- ~JJ v- V ct^+VSF = Vf + VSF, (17) 
where the force density / contains the advective, viscous and stochastic contributions to the fluid 
dynamics. An important feature of this formulation is that the density p in the usual Navier-Stokes 
equation is now replaced by the effective density operator 

Pca = P I + m e VSJV, (18) 
where J is the identity operator or matrix. Notice that the effective density operator for incom- 



pressible flow is not pi + m e SJ as one might naively expect based on ( 11 ). The distinction between 
SJ and VSJV is important, and leads to a well-known but surprising difference between the short- 
time motion of a particle immersed in a compressible versus an incompressible fluid [50] . When 
there are many particles present, the effective inertia tensor is generalized straightforwardly by 
summing the added inertia over the particles, 

p eS = P I + Y,{m e ) i VS i J i V. (19) 



Equation (17) together with the no-slip condition q = Jv gives a closed set of equations for 
v and q without any constraints. We can use this unconstrained formulation to simplify anal- 
ysis of the properties of the coupled fluid-particle problem. In particular, in Appendix [Bj the 
constraint-free form is used for showing fluctuation-dissipation balance in the stochastic setting. 



For numerical purposes, however, (17) may be difficult to use directly because of the appearance of 
the cumbersome operator VSJV. However, as we explain in Appendix [AJ with periodic boundary 
conditions it is possible to efficiently invert the operator p cS using Fourier transforms and thus 
obtain a closed-form equation for d t v suitable for numerical implementations. This relies on the 
fact that for a large d-dimensional periodic system JVS is a constant multiple of the d x d identity 
matrix. 

E. Fluctuation-Dissipation Balance 

So far, we considered the equations of motion for the fluid-particle system ignoring thermal 
fluctuations. In Appendix [B] we formally demonstrate that in order to account for thermal fluctu- 
ations in a manner that preserves fluctuation-dissipation balance it is sufficient to add the usual 
Landau-Lifshitz stochastic stress (k B Tr/) 1/2 (W + W T ) to the viscous stress tensor in <x, without 
adding any stochastic forces on the particle. The key physical insight is that the fluid-particle 



coupling is non-dissipative, as demonstrated in Section II C and the only dissipation comes from 
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the viscous terms. 

Fluctuation-dissipation balance here means that at thermodynamic equilibrium the particle- 
fluid system is ergodic and time-reversible with respect to the Gibbs-Boltzmann distribution 



Z~ x exp (—H/k B T), where the "Hamiltonian" H given in (16) is to be interpreted as a coarse-grained 

free energy. Since u = Jv is not an independent degree of freedom, we can formally write the 

Hamiltonian in terms of the degrees of freedom of the system as a sum of potential and kinetic 

energy. The total kinetic energy includes, in addition to the kinetic energy of the fluid J (p/2) v 2 dr, 

a kinetic energy contribution due to the motion of the particle, 

E„ = m e — = m„ = — - / v ■ (Su) dr — — - / v T SJv dr, 

p 2 2 2 J K ' 2 J 

where use was made of the adjoint condition ([2]). This leads to the coarse-grained Hamiltonian 

/ii T <? Tn f v 2 1 f 

-^dr + p J -dr + U(q) = -j v T p cH v dr + U (q) . (20) 

Note that it is not necessary here to include an entropic contribution to the coarse-grained free 
energy because our formulation is isothermal, and we assume that the particles do not have internal 
structure. 



We emphasize that the form of H (q, v) in (20) is postulated based on physical reasoning rather 
than derived from a more refined model. Atzberger |27j provides a careful and detailed discussion 
of how one might eliminate the particle velocity u by performing an adiabatic elimination, starting 
from a frictional coupling model in which the particle velocity is independent from the fluid velocity. 
The starting frictional coupling model, however, as we discussed earlier, involves an arbitrary 
frictional force parameter that is yet to be given a microscopic interpretation. We believe that the 
consistency of our inertial coupling model with general thermodynamic principles and deterministic 
hydrodynamics is sufficient to adopt the inertial coupling model as a consistent coarse-grained 
model without having to justify it from "first principles" |51j . 



The fact that fluid-particle coupling conserves the Hamiltonian (see Section II C ) and is therefore 
non-dissipative is a crucial component of fluctuation-dissipation balance. However, this is not 
sufficient on its own. An important additional requirement is that the phase space dynamics should 
be incompressible, which means that the dynamics preserves not just phase-space functions of H, 
but also preserves phase-space measures such as the Gibbs-Boltzmann distribution. As discussed 
in Atzberger [27] , even for the case of a neutrally-buoyant particle an additional "Ito" or "thermal" 
drift term needs to be added to the velocity equation to ensure fluctuation-dissipation balance. 
This term has the form of an additional contribution to the stress tensor 

<r th = [S (k B T)} I. 



17 



The physical origin of this term is the Kirkwood kinetic stress due to the thermal motion of the 
particle lost when eliminating u as a degree of freedom. Another way to interpret this term is that 
it adds the contribution of S (k B T) to the pressure tt coming from the particle. For incompressible 
flow, this simply changes the pressure but does not change the dynamics of the velocity field since 
the projection "P eliminates the scalar gradient term V • <r th = VS (k B T). In Appendix [B] we argue 
that for non-neutrally buoyant particles there is also no need to include an additional thermal drift 
term for periodic boundary conditions. 



1. Equipartition of Energy 

The fact that the Gibbs-Boltzmann distribution is separable in v and q and that the Hamiltonian 
is quadratic in v means that the fluctuations of velocity are Gaussian with covariance (vv*) = 
(k B T) p~g. The fluctuations of the particle velocity have variance 

(u 2 ) = Trace [J (vv*) S] = (k B T) Trace [Jp^S] 
For a single particle immersed in a periodic incompressible fluid in d dimensions, using the relations 



(|A5|) and (|A3|) derived in the Appendix, we can simplify 

) 



k B T 

Irace 



,k B T 



= d^, (21) 



P 

where AV = dAV/ (d — 1) and rh — m e + dm f / (d — 1). 

The result (u 2 ) — d (k B T) /rh should be compared to the corresponding result for a compressible 
fluid |38j, (u 2 ) = d(k B T) /m, which follows from the usual equipartition principle of statistical 
mechanics. When incompressibility is accounted for, a fraction of the equilibrium kinetic energy 
is carried in the unresolved sound waves, and therefore the apparent mass of the particle is rh 



and not m = m e + m f , as we verify numerically in Section IV C 1 It is reassuring that our model 
equations reproduce the result for rigid particles. This suggests that the model introduced here 
can be used to study more complicated questions such as the effect of multi-particle interactions 
on (u 2 ) in semi-dilute to dense colloidal suspensions 



2. Brownian Dynamics Limit 

In the long-time limit, the motion of a single free particle immersed in a fluid looks like simple 
Brownian motion with a diffusion coefficient that can be defined from the long-time mean square 
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displacement of the particle. In general, one defines the diffusion tensor 

v , 1 /[g(t)-g(0)] [q(*)-q(0)] T \ 

While we will not demonstrate this here, it can be shown that, for very small Reynolds numbers, 
neither sound effects (compressibility) nor inertial effects affect the diffusion coefficient of a single 
particle [S3] . 

For a single particle, x«, can be obtained from the mobility tensor fi via the Stokes-Einstein 
relation Xoo — (k B T) fi. Here the mobility operator is defined from the relation u = (J,F , where 
F is a constant force applied on the particle and u is the steady-state velocity of the particle 
under the effect of this applied force Specifically, let v = -C~ 1 SF be the solution to the 
steady-state Stokes problem, 



i]V 2 v = - Vtt + SF , such that V • v = 0, 

with appropriate boundary conditions. Here £ _1 is simply a short-hand notation for the linear 
Stokes solution operator. The velocity of the particle is then 

Mo = Jv = - (JC~ X S) F = fj,F . 
For a periodic system, the mobility does not depend on the position of the particle, fi(q) = fi 
by translational symmetry, and one can show using Fourier transform techniques that the different 
dimensions are not coupled. In three dimensions, this defines an effective hydrodynamic radius R H 
for the particle 

T 1 T a T 

&TrrjR H na 

where a is the physical size of the particle defined through the extend of the kernel function 
5 a (Ar), and a is some constant that depends on the particular choice of the kernel function and 
also the system size. It is important to note here that the diffusive motion of the particle is 
entirely determined by its coupling to the fluid and is not an input parameter to the method. This 
reflects the physical relationship between fluid velocity fluctuations (viscosity and temperature) and 
diffusion coefficient, as encoded in the Stokes-Einstein relation. If it is desired to tune the diffusion 
coefficient of the particles to match some, for example, experimental value, one may increase it by 
adding a random slip as described in Appendix [Bj 

For non-periodic boundaries the mobility in general depends on the position of the particle 
relative to the boundaries. For multi-particle suspensions the multi-dimensional mobility tensor 
M (Q) = {fJ-ij} depends on the positions of all particles Q — {q t }, and the block corresponding to 
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particle pair i and j is 

In this case, the precise statement about the Brownian dynamics of the particles is that in the 
overdamped limit, in which momentum diffuses much faster than the particles, the motion of the 
particles at the diffusive time scale can be described by the fluid-free Ito stochastic differential 
equation [27] 

Q = MF + ^2k B T M^W + k B T \Jq- m J , ( 22 ) 

where W is a collection of independent white-noise processes and F (Q) are the forces applied on 
the particles. 



The validity of this very common Brownian dynamics approximation (22) rests on the separation 
of time scales between the diffusion of momentum and the motion of the particles. This is usually 
estimated by the Schmidt number defined as 



v v arf 



Xoo (k B T) n apk B T 

The Brownian dynamics limit may be used only if S c ^> 1, which is typically the case for micron- 
scale particles but not for nano-scale particles. However, it is important to understand that the 
Brownian dynamics limit only strictly applies to very long times, once momentum has diffused 
the length of the system L, and not just the typical size of a single particle or even the typical 
separation between nearby particles. 

The short-time motion of particles immersed in the fluid is known to be very strongly affected 
by momentum diffusion and by inertial effects [53] . Since our method does not make the Brownian 
dynamics approximation, it is able to produce both the correct short-time and long time features 
of the Brownian motion of a particle, where "short" refers to time scales after sound waves have 
decayed. This observation was made for the Stochastic Immersed Boundary Method in Refs. 
[544 [55] , and here we demonstrate that our Inertial Coupling Method also correctly captures the 
known physical effects of particle inertia. 

III. SPATIO-TEMPORAL DISCRETIZATION 

In this section we describe our second-order spatio-temporal discretization of the equations 
of motion ( 7|8|9| ). Our spatio-temporal discretization is based on the deterministic Immersed 



Boundary Method (IBM), and in particular, on the deterministic second-order temporal integrator 
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presented in Ref. [56]. For the fluctuating fluid solver, we use the second-order staggered-grid 
spatial discretization of the fluctuating Navier-Stokes equations described in detail in Ref. |39j . A 
notable feature of the fluid solver we employ is that it handles the viscous terms semi-implicitly and 
is stable for large time steps. Furthermore, for the fluctuating Stokes equations, our fluid solver 
produces the correct spectrum of the velocity fluctuations for any time step size |39|. 

There are two key novel features in our incompressible inertial coupling algorithm from those 
previously developed. Firstly, our algorithm includes the effects of particle excess inertia in a man- 
ner that strictly conserves momentum and is second-order deterministically for smooth problems. 
Secondly, we focus our initial development on systems with periodic boundaries only, allowing the 
use of the Fast Fourier Transform (FFT) as a linear solver for the time-dependent Stokes equa- 
tions. This greatly simplifies the implementation of the algorithm and allows us to use Graphics 
Processing Units (GPUs) for very efficient parallelization of the algorithm. 

For neutrally-buoyant particles our temporal discretization is exactly that described in Ref. 
|56| with the fluid solver replaced by that described in Ref. |39j . This simplified algorithm is 
already implemented in the IB AMR software framework [57] , an open-source library for developing 
fluid-structure interaction models that use the immersed boundary method. Note that IBAMR 
can handle non-periodic boundary conditions using a preconditioned iterative solver for the time- 
dependent Stokes equations [58] . In this paper we focus on describing the additional steps required 
to handle excess inertia and to use FFTs as a linear solver. 

The majority of our presentation will focus on a single particle coupled to a fluctuating fluid. 
Only small changes are required to handle multiple particles by simply summing the single-particle 
term over the different particles. As we explain in more detail shortly, the error introduced by 
superposing the single particle solutions to solve the multi-particle system is small if the kernels of 
the different particles are not overlapping, which in practice means that there are at least 3 grid 
cells between the centroids of the particles. 



A. Spatial Discretization 



Our second-order spatial discretization of the equations ( 7|8|9 ) is based on standard techniques 



for incompressible flow and the immersed boundary method [30] , as described in more detail in, for 
example, Ref. [56]. In the spatially-discretized equations, the same equations as for the continuum 
apply, but with the interpretation that the velocity v is not a (random) vector field but rather a 
finite- volume discretization of that field [59, 60J. We use a uniform Cartesian staggered-grid spatial 
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discretization of the incompressible Navier-Stokes equations, as described in more detail in Ref. [39J . 
In the staggered discretization the control volume grid associated to each component of velocity 
is shifted by half a grid spacing along the corresponding dimension relative to the pressure grid. 
In the discrete setting, the various continuum operators acting on vector fields become matrices. 
The spatial discretization of the differential operators, notably the discrete gradient, divergence 
and Laplacian operators, is described in detail in Ref. [59] . 

Application of the local averaging operator J, which is a convolution operator in the continuum 
setting, becomes a discrete summation over the grid points that are near the particle, 

Jv= ^2 <$><>, (q-r k )v k , 

fceEgricl 

where r k denotes the center of the control volume with which v k is associated, and 4> a is a function 
that takes the role of the kernel function S a . We follow the traditional choice [30J and do the local 
averaging independently along each dimension, 

d 

4>a (<? r k ) = Y[ 4>a [q a - (r k )J , 

a = l 

which improves the isotropy of the spatial discretization (but note that the local averaging is not 
rotationally invariant). As a matrix, the local spreading operator S = (AVfY 1 J* is a weighted 
transpose of J, 

(SF) k = (AV f y 1 cj )a (q-r k )F, 

where AV f = AxAyAz is the volume of the hydrodynamic cell 1 . 

The discrete kernel function 4> a was constructed by Peskin [30] to yield translationally-invariant 
zeroth- and first-order moment conditions, along with a quadratic condition, 

fcggrid 

£ 0?~ r k )4> a (q- r k ) = 

fcGgrid 

£ 4>t(q-r k ) = AV^ 1 = const., (23) 

fcggrid 

independent of the position of the particle q relative to the underlying (fixed) velocity grid. En- 
suring these properties requires relating the support of the kernel function to the grid spacing, 
that is, making a ~ Ax. This means that the size and shape of the particles is directly tied to 
the discretization of the fluid equations, and the two cannot be varied independently, for example, 
simulating the motion of a "spherical" particle requires choosing the same grid spacing along each 

1 The cell volume AVf is introduced here because the fluid kinetic energy appearing in the discrete Hamiltonian is 
AVf X^feggrid P v 1/2 and therefore dH/dv is not the functional derivative pv as in the continuum (see Appendix 
IB]) but rather the partial derivative AVf pv. 
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dimension, Aa; = Ay = Az. This is a shortcoming of our method, but, at the same time, it is 
physically unrealistic to resolve the fluid flow and, in particular, the fluctuations in fluid velocity, 
with different levels of resolution for different particles or dimensions. 



The last condition (23), imposed for unrelated reasons by Peskin [30], is especially important 
in our context since it implies that the particle volume AV = (JS)' 1 will remain constant and 
independent of the position of the blob relative to the underlying grid. The function with minimal 



support that satisfies (23) is uniquely determined [30J. In our numerical experiments we employ this 
three-point discrete kernel function, which means that the support of </> extends to only three grid 
points along each dimension (i.e., 3 d discrete velocities are involved in the averaging and spreading 
operations in d dimensions) |38j . This particular choice gives AV = 8AVf in three dimensions and 
Ay = 4AVf in two dimensions. The narrow kernel improves the computational efficiency on the 
bandwidth- limited GPU, as detailed in Ref. |38j. 

1. Translational Invariance 

In the continuum derivation, obtaining a closed-form pressure-free velocity equation relied sen- 
sitively on the fact that for a large continuum system 

JVS= ^— -AV-'I, (24) 
a 

see Eq. (A3) in Appendix [Aj Ideally, we would like the spatial discretization to have the additional 



property that, for periodic boundary conditions, JVS be invariant under translations of the particle 
relative to the underlying fluid grid. This is not ensured by the Peskin operators, which are 
constructed without any reference to the fluid equations and the particular form of the discrete 
projection or the discrete viscous dissipation. In fact, in the traditional immersed boundary method 
|30| . a centered discretization of the velocity was used, which implies a very different form for the 
discrete projection operator V and required the introduction of an additional "odd-even" moment 
condition not strictly necessary with a staggered discretization. 

Numerical experiments suggest that for the staggered grid discretization that we employ, the 



continuum identity (24) is obeyed to within a maximal deviation of a few percent, 

aF' (q) = JVS « AV' 1 1 (25) 
for a sufficiently large periodic system, where AV = dAV/ (d — 1) is a modified volume of the blob. 
For a periodic three-dimensional system of N x x N y x N z hydrodynamic cells, we expect to see 



deviations from ( 25 ) if one of the grid dimensions becomes of the order of the kernel width (which 



is 3 cells for our spatial discretization). This is because a particle then becomes affected by its 



23 



nearby periodic images. For a grid size N x x N y x 1 cells we expect to obtain two-dimensional 



behavior [see Eq. (A4)]. 

In the left panel of Fig. [T] we show numerical results for the average and maximum deviation of 
AV (JVS) from the dx d identity matrix, 

5I(q) = -^-jAV(JVS)-I, 
as we vary the position of the particle q relative to a cubic Eulerian grid of size N x = N y — N z = N 
cells. Specifically, we show the average diagonal value, the maximum diagonal element, and the 
maximum off-diagonal element of SI (q). For all but the smallest systems the diagonal elements are 
smaller than 5% and the off-diagonal elements are on the order of 0.1%. For smaller system sizes 
there are visible finite-size effects due to interactions with periodic images. For comparison, we 
also show the corresponding two-dimensional results for a square grid of N x — N x = N cells. The 
finite-size effects are more pronounced in two dimensions due to the slower decay of the Green's 



function for the Poisson equation, but for systems larger than N = 32 cells we find ( 24 ) to hold to a 
percent or so. In the left panel of Fig. [T]we also show the average diagonal values (SI xx +SI yy )/2 and 
8I ZZ for non-cubic systems, illustrating the change from three-dimensional to the two-dimensional 
behavior as N e — > 1. 

Given these numerical observations, we will make an approximation and assume that JVS is 



a constant multiple of the dx d identity matrix independent of q, and (25) holds as an equality. 
While this is, in principle, an uncontrolled approximation, we will correct for it in order to en- 
sure strict momentum conservation and thus strict adherence to fundamental physical laws. This 



approximation allows us to write the discrete equivalent of ( A5 ) 

'SJJ 



(pl + meVSJVy 1 w p- 1 ll-^^VSJV } . (20) 



When there are multiple particles in the system we simply sum the second term on the right hand 
side over all particles, which is an approximation even in the continuum setting. This approximation 
relies on the assumption that J{PSj ~ for two particles i and j that are far away from each other. 
To test the validity of this approximation, in the right panel of Fig. [T] we show the maximum 
diagonal and off-diagonal value for J{PSj for two particles a distance I away from each other, for 
both q i — q j = I (1,0,0) and q i — q 3 = I (\/3, y/3, y/S) /3, in three dimensions. If I = then J {PSj 
approaches J/3, so we normalize the value of J {PSj by a factor of SAV/2. We see that for 

I > 5Ax, the maximum normalized value of J{PSj is less than 0.01. As seen in the figure, JiSj 
vanishes identically if the kernels of particles i and j are disjoint, which is always the case for 
I > 3Ax. 
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Figure 1: (Left panel) Translational invariance of the approximation dAV (J'PS) / (d — 1) « I over a set of 
10 3 positions of the particle relative to the underlying fluid grid. For a periodic three-dimensional system 
of iV 3 cells, we show the average diagonal value of 81 (black squares), the maximum diagonal element 
of SI (magenta stars), as well as the maximum off-diagonal element of 51 (cyan pluses). We also show 
the average diagonal value of 51 for a two-dimensional system of iV 2 cells (red circles), as well as the 
average value (51 xx + 5I yy )/2 (green diamonds) and 51 zz (blue triangles) for a three-dimensional system 
of 32 x 32 x N cells. (Right panel) Maximum diagonal and off-diagonal value of 3AT^ (JiPSj) /2 and the 
maximum diagonal value of AV (JiSj) for two particles i and j a distance I apart along the (1,0,0) or 
(1, 1, 1) direction, in a three-dimensional periodic box of 32 3 cells. 

Another operator that is important for the long-time diffusive (Brownian) motion of the particle 
and enters in the overdamped equations, as explained in Section |IIE2[ is the discrete mobility 
operator for a single particle, 



C19J 



where C 1 denotes the discrete Stokes solution operator. For a single particle in a periodic domain, 

oximately t 

m(<z) ~ 



we numerically find that (q) is approximately translationally-invariant to within a few percent 

1 



6nr]R H ' 

where the effective hydrodynamic radius for an infinite system is numerically extrapolated to be 
R H « 0.91Aa; for a uniform grid spacing Ax and the three-point kernel [38]. This is consistent with 
previous results for a cell-centered discretization of the Navier-Stokes equations |61| . By using the 
Peskin four-point kernel [30j instead of the three-point discrete kernel function the translational in- 
variance of the spatial discretization can be improved, however, at a potentially significant increase 
in computational cost. A more systematic investigation of different choices for the discrete kernel 
functions has been performed by Mori however, these types of investigations have yet to be 
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carried out within the "blob" particle approach. We have found the cheap three-point function to 
perform quite well in our tests and use it exclusively in this work. 

For non-periodic systems, one must generalize the definition of J and S in the case when the 
particle overlaps a physical boundary [63J. Even if a particle does not overlap a boundary, however, 
it will feel the boundary hydrodynamically and therefore both JVS and JC~ X S will depend on the 
proximity of the particle to physical boundaries. Implementing our algorithm in such cases may 
require first pre-tabulating the values of these d x d matrices for different positions of the particle 
relative to the boundaries. 

B. Temporal Discretization 

In this section, we describe how to integrate the spatially-discretized equations in time and 
update the fluid and particle velocities and particle position from time nAt to time (n + 1) At, 
where At is the time step size, which can in principle be adjusted dynamically but we will assume 
it is kept fixed. We will use a superscript to denote the time level at which a given quantity is 
evaluated, for example, q n+ i will denote a mid-point estimate for the position of the particle at 
time (n + |) At. Similarly, F n+ ^ = F (q n+ ^ will denote the force on the particle (due to external 
sources or other particles) evaluated at the position q n+ h and, implicitly, at time (n + |) At in the 
case of a time dependent force. 

For a neutrally-buoyant particle, m e = 0, the fluid momentum equation is coupled to the particle 
position only through the forcing term SF. In the deterministic setting, Griffith and Luo [56j have 
developed a second-order splitting scheme for integrating the spatially-discretized equations in 
time. The fluid solver in this scheme is very similar to the predictor-corrector scheme employed in 
the stochastic setting in Ref. [39]. The temporal discretization that we present next is based on 
replacing the fluid solver in Ref. |56j with that in Ref. [39] . at least for the simpler case m e = 0. 
The main difficulty is including the additional inertia from the particles in the fluid momentum 
update in a computationally-efficient manner. 

In Ref. [38j a first-order splitting algorithm was developed for the case of a compressible 
fluid. This type of algorithm is similar to the original projection algorithm of Chorin [63] for 
incompressible flow and can be summarized as follows. Update the fluid first without accounting 
for the force A exerted by the particle. Then, solve for the value of A that, when applied as 
a correction to the fluid update, exactly imposes the no-slip condition. Extending this type of 
approach to be higher than first order accurate is known to be difficult from the literature on 
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incompressible flow [65], due to the fact that the splitting introduces a commutator error. Here we 
follow a different, though related approach, which allows to construct a more accurate algorithm 
for viscous-dominated flows. 

Our temporal scheme will be based on the following approach, which can be shown to be 
second-order deterministically by a Taylor series expansion of the temporal local truncation-error: 

1. Estimate the position of the particle at the midpoint to leading order, 

q n+l ^ g n + ^*J»„». (27) 



m e SJ ( v ■ -^-J ) v 
oq 



2. Update the fluid velocity based on (11) using a second-order algorithm, while keeping the 
particle positions fixed at the midpoint estimates, 

(pi + m e S n+ ? J"+*) — + Vtt"+? = -V- (pvv T + cr) n+h +S" +1 -F n+ K 

L 

(28) 

subject to V • v n+1 — 0. Here a second order Runge-Kutta [42J or Adams-Bashforth [56] 
scheme can be used to evaluate the fluid momentum fluxes to at least second-order accuracy, 
denoted generically here by superscript n + | . 

3. Update the particle position using a second-order midpoint estimate of the velocity, 

g n + l = Q n + At jn+1 ( v n + l + _ (29) 

Observe that the above scheme never actually uses the particle velocity u, although one can and 
should keep track of the particle excess momentum m e u and update it whenever the fluid momen- 
tum is updated, to ensure strict conservation of momentum. Also observe that during the fluid 
update we fix the particle at its midpoint position q n+ ?. 



1. Velocity Update 

The most difficult step in the time stepping algorithm summarized above is the momentum 
(velocity) update, step[2j In order to update the velocity of the fluid we need to calculate the fluid 
momentum change due to viscosity and thermal fluctuations and also the momentum exchange 
with the particle, all to second order in time. Our scheme is based on solving for the values of the 
Lagrange multipliers 7r n+ 2 and X n+ ^ such that at the end of the time step both the incompressibility 
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and the no-slip constraints are satisfied, 

v n + l _ n 



At 



-V- (pvv T + <r) n+ i -S n+ i\ n+ * 
m e u n + AtF n+ s + At A n+ 3 




(30) 



A correction A«" + j is included to account for the fact that the no-slip condition is not correctly 
centered since J and v are evaluated at different points in time, as we explain shortly. 

If we multiply the particle velocity update by S n+ ^ and add it to the fluid equation, and use 
the no-slip constraint including u n = J n ~^v n + Au n ~i , we get 

; pi + m e S n+ ^J n+ ^ A + Vtt"+5 = -V • (pvv T + <r)" +5 + S n+l -F n+1 - 



A/ 



At 



(31) 



which is to be compared to the desired (28). The penultimate term in the above equation can be 



seen as a discretization of the kinetic term 
This is consistent with the continuum equations but it yields only first-order and not second-order 



u ■ J ! c 



accuracy in At because it is not centered at time level n + \ . This reduction of the accuracy comes 



because the no-slip constraint in (30) uses the midpoint instead of the endpoint position of the 
particle. 

The above discussion shows that setting Au n+ i = results in a first-order scheme if m e ^ 0. To 
get second-order accuracy, we need to apply a nonzero correction to the no-slip constraint. Imposing 
the no-slip constraint at the end of the time step, u n+1 = J n+1 v n+1 , leads to a formulation that is 
implicit in both q" +1 and v n+1 , which is difficult to implement in practice. Instead, we can center 
the no-slip constraint as 

' s (v n+1 +v n ) = -(«." 1 + J>" i = - ( J" 1 ' 1 + At/" 1 i - J'V 



2 v ' 2 v 

which gives the no-slip centering correction 

Au"+3 = (j"+5 - J") u n . (32) 

This correction for the no-slip constraint is simple to implement with only one additional local 
averaging operation to evaluate J n v n . Note that we purposely used J n v n instead of u n here since 
in our formulation, and also in our algorithm, u n+1 is only used as an intermediate variable. A 



Taylor series analysis shows that using (32) makes the last line in (31) a centered second-order 
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approximation of the kinetic term 



Am"- 2 



At 



At 



u ■ — J ) u 



0(Ai 3 



A Taylor series analysis confirms that using the no-slip correction (32) leads to a second-order 
algorithm for updating the position of the particle and the velocity of the fluid. 

In order to avoid one more additional local averaging operation (which requires an irregular 
memory access pattern and is thus costly, especially in a GPU-based implementation) we can set 
Au n+ ? = 0. We are primarily concerned with viscous-dominated (low Reynolds number) flows, for 
which the kinetic term (u ■ dJ/dq) v is small (quadratic in v, just like the advective term »-V«) and 
can be approximated to first order without a significant reduction in the overall accuracy of the 
method. As we explain in Appendix [cj in our scheme Au n+ ? contains an additional higher-order 
correction (0(Ai 3 ) for smooth flows) that arises solely due to the implicit handling of viscosity. 



2. Semi-Implicit Discretization of Viscous Terms 

During the fluid update the particle position remains fixed at q n+ ^. For notational simplicity, 
in the remainder of this paper we will sometimes drop the time step index from J and S; unless 
otherwise indicated, they are always evaluated at q n+ i . 



Following Ref. [39], our second-order implementation of the velocity update (31) treats the 
viscous term semi-implicitly and the remaining terms explicitly, 

-V • (pvv T + a) n+h = |V 2 (v n+1 + v n ) + V • £ n - V • (pvv T ) n+h . 



The spatial discretization of the stochastic flux is [39] 

*- = {-^)" lw " + iw " f] - 

where AV f = AxAyAz is the volume of the hydrodynamic cells, and W" is a collection of i.i.d. 
unit normal variates, generated independently at each time step on the faces of the staggered 
momentum grid. To approximate the advective fluxes to second order in At, one can use either 
the predictor-corrector method described in Ref. [39J or, more efficiently, one can use the classical 
(time lagged) Adams-Bashforth method [56] 

V • (pvv T ) n+h - • (pvv T ) n - • (pvv^ 1 - 1 . (33) 
For viscous-dominated (small Reynolds number) flows, one can also approximate the advective 
terms to first order only without a significant reduction of the overall accuracy for reasonably large 
time steps. 



Referring back to Eq. (28), we see that updating the fluid momentum semi-implicitly requires 
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solving the linear system 
' pI + m e S n+ ^J n+ * 



At 



-V • (pvv 



T\»+5 



V ■ 17 



pI + m e S n+ 2j 



Af 



S n+ 2 p n +o 



m e SJ I f • — J ) v 



(34) 



If m e = 0, we can solve the linear system (34) for the unknowns v n+1 and 7r n+ 2 using a preconditioned 



iterative solver [58], as explained in more detail in Ref. [39]. For periodic systems the system (34) 



can be solved easily by using a projection method together with FFT-based velocity and pressure 



linear solvers. For non-neutrally-buoyant particles, however, solving (34) requires developing a 



specialized preconditioned Krylov method. Here we develop an approximate solver for (30) that 



only requires a few FFTs, and will be shown in Section IV A to give nearly second-order accuracy 
for a wide range of relevant time step sizes. 

Our approach consists of splitting the velocity solver into two steps. In the first step, we ignore 



the inertia of the particle, i.e., delete the m e SJ term in (34), and solve for a provisional velocity 
i>" +1 and pressure 7f n+ 2, 

( Xt 1 ~ I v ') + v * n+h = (IV + 2 v v * + v ' s " + spn+h ~ v ' O^T^ . ( 35 ) 

subject to V • v n+1 = 0. If m e — 0, this completes the fluid solve and setting v" +1 = v n+1 gives 



us second-order accuracy for the viscous and stochastic terms |42j . If m e ^ 0, we need to find a 
velocity A»" + s = v n+1 ~v n+1 and pressure correction A-7r n+ 5 = n n+ i — 7r n+ 5 that takes into account 



the inertia of the particle. We do this by splitting the linear system (30) into two equations, (35) 
for the unperturbed velocity field, and 



V ( Att" 



-u 



n + l 



+ F n 



At At 

u n+1 = j(v n+1 + Av r ' 



Am" 



(36) 
(37) 
(38) 



V • (Au™+2) = 0, 

for the perturbed field. This gives a linear system of equations for the unknowns v n+1 , u" +1 , Av n+ i , 
jr"+h , A-K n+ i , and A n+ s . We explain how we solve this linear system of equations in Appendix [c] for 
periodic boundaries using Fourier Transform techniques. Here we simply summarize the resulting 
algorithm, as implemented in our code. In Appendix [D] we give a summary of a similar algorithm 
for compressible flow, which our code also implements. 



C. Summary of Algorithm 

1. Estimate the position of the particle at the midpoint, 

and evaluate the external or interparticle forces F n+ ^ . 

2. Solve the unperturbed fluid equation 

i> n+1 + i, n \ ■+■ V . I 



• (pvv T ) n - -V • ( P ™ T ) r 



5™+ 2 _P™+2 ; 



V-v n+1 = 0, 

using a projection algorithm and FFTs to diagonalize the Laplacian operator. 

3. If m e = 0, set v n+1 = v n+1 and skip to step|9| 

4. Evaluate the slip correction 

5u"+? = - J") v" + ^J"-5V 2 (a« b -*) 

and the change of the particle excess momentum 

Ap = m e (u n - J n+ ^i> n+1 - Su n+ ^ . 

5. Calculate the fluid velocity perturbation due to the excess inertia of the particle 

777 t 

A* - -r^r ;VSAp, 

pyrrif + m e ) 

using FFTs to implement the discrete projection "P, where m f = dpAV/ (d— 1). 

6. Account for the viscous contribution to the velocity perturbation by solving the syste 
pi- ^??V 2 J Av n+ ? + AtV (A7r n+ *) = S n+ ^Ap-m e J n+ ^Av 



V • ( Av n+ * J = 

using a projection algorithm and FFTs to diagonalize the Laplacian operator. 

7. Update the fluid velocity 

v n + l = {3" + 1 + Av"+5. 

8. Update the particle velocity in a momentum-conserving manner, 

u n + l = jn+I /-„+! + A ~N + ^n+i _ 
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9. Update the particle position, 



qn+ i = q n + + „») . (48) 



We note that the full slip correction (43) is only required if m e /m/ is large and the Reynolds 



number is large. For sufficiently small Reynolds numbers (viscous-dominated flows) we can neglect 
the quadratic advective term and only keep the linear term, and set 

5u ™+h = ^j-^v 2 (A^-s) . (49) 
We can also set Su n+ i = 0, and obtain a first-order algorithm that does not require any time lagging 



and has improved stability for very large time step sizes. We compare the three options (43), (49), 



and 5u n+ a = numerically in Section IV A The remainder of the algorithm is not affected by the 
choice of the slip correction Su n+ K 



IV. RESULTS 



In this section we validate and test the performance of the algorithm summarized in Section 



III C With periodic boundary conditions the velocity and the pressure linear systems in the incom- 
pressible formulation decouple and Fast Fourier Transforms can be used to solve them efficiently 
|30| . More precisely, a splitting or projection algorithm can be used to solve the Stokes (velocity- 



pressure) system (40). We first solve the velocity equation (40) without the gradient of pressure 
term (this is a Helmholtz equation) using a Fourier transform to diagonalize the discrete Lapla- 
cian. Then, we project the solution onto the space of divergence free vector fields by subtracting 
a pressure gradient term. The pressure is a solution of a discrete Poisson equation, which can also 
efficiently be computed using Fourier transforms. Note that it is possible to generalize our algo- 
rithm to non-periodic systems by using the fluid solver developed by one of us [58] and employed 
in Ref. |39| . at least for the case of neutrally buoyant particles, m e = 0. For m e ^ new iterative 
solvers for the Stokes subproblem need to be developed. 

We have parallelized the algorithm to run efficiently on Graphics Processing Units (GPUs), as 
explained in more detail in Ref. [35] • Our public domain implementation [3H] is written in the 
CUD A programming environment, and is three-dimensional with the special case of N z = 1 cell 
along the z axes corresponding to a quasi two-dimensional system. In our implementation we create 
one thread per cell, and each thread only writes to the memory address associated with its cell and 
only accesses the memory associated with its own and neighboring cells. This avoids concurrent 
writes and costly synchronizations between threads, facilitating efficient execution on the GPU. 
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For incompressible flow, our present GPU implementation is specific to periodic systems, and uses 
the NVIDIA FFT library as a Poisson/Helmholtz solver. 

The stability and accuracy of our spatio-temporal discretization is controlled by the dimension- 
less advective and viscous CFL numbers 

a= AF- (50) 

where V is a typical advection speed, which may be dominated by the thermal velocity fluctuations 
or by a deterministic background flow. Here we always use the same grid spacing along all dimen- 
sions, Ax — Ay = Az. The strength of advection relative to dissipation is measured by the cell 
Reynolds number r = a/0 — VAx/u. Note that for compressible flow (see Ref. [39] and Appendix 
[P] ) there is a sonic CFL number a s — cAt/Ax, where c is the speed of sound. 

The explicit handling of the advective terms places a stability condition a < 1, in fact, for 
a > 1 a particle can move more than a hydrodynamic cell during a single time step and this causes 
not only stability but also implementation difficulties. It is not hard to see that in the absence 
of advection our semi-implicit discretization of viscosity is stable for any value of /3, however, it 
is only by keeping (3 < 1 that we can ensure the dynamics of all or at least most fluid modes is 
resolved [42J. We consider a temporal integrator to be "good" if it produces reasonably-accurate 
results with a time step for which at least one of a or /3 is close to 1/2. Typically, flows at small 
scales are viscous dominated (r <C 1) so that the time step is primarily limited by (3 and not by a. 

A. Deterministic convergence tests 

Our first test is a check of the deterministic order of accuracy of the temporal integrator. 
Based on local truncation error analysis we expect that the temporal integrator summarized in 



Section IIIC is formally second-order accurate. For small Reynolds numbers, we expect to see 



nearly second-order accuracy in practice even if we use the slip correction (49) instead of (43). 
We also recall that for m e ^ we made an uncontrolled approximation in assuming that JVS is 
translationally-invariant, which is only accurate to about a percent for the three-point Peskin local 
averaging and spreading operators. This approximation leads to another error in imposing the 
no-slip condition, which we expect to lead to first-order accuracy for very small time step sizes. 

As a test of the temporal accuracy, we study the deterministic motion of a particle in an 
centrally-symmetric harmonic potential V(r) — kr 2 /2, where r is the distance from the origin and 
A: is a spring constant. In these tests we keep the spatial discretization (and thus the blob particle 
shape) fixed and only change the time step size At. We start the particle from rest at a certain 
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Figure 2: Average error (51) over a short deterministic trajectory of a particle initially at rest and subse- 



quently moving under the action of a harmonic potential. (Left panel) Small Reynolds number, Re « 0.02, 



with no-slip correction (49 ). Several values of the excess particle mass m e relative to the mass of the dragged 
fluid rrif (symbols, see legend) are shown, including neutrally-buoyant particles (m e = 0). Expected error 
decay rates for a first and second-order scheme are shown with lines. (Right panel) Comparison of the three 



choices for the slip correction, (43 1 (viscous and centering corrections), (49) (only viscous correction), and 
(43) (no correction), for Reynolds numbers Re sa 2 (full symbols) and Re w 20 (shaded symbols). The excess 
mass is m e = 10m/. 



distance r from the origin and then release it. The particle will perform damped oscillations under 
the influence of the spring and viscous friction. We look at the error in the position of the particle 
q(t) defined as the average of the difference between the position for time steps At and At/ 2 over 



a certain number of time steps N s , from the initial time to a time T = N s At, 



E(At) 



2n- 



At 



(51) 



q A t (nAt) - q At/2 

For a numerical scheme with order of accuracy p this error should behave as E — O (At p ) for 
sufficiently small At. 

We perform this test for several choices of the density of the particle relative to the fluid, as 
measured via the ratio m e /mf, including particles less dense than the surrounding fluid (negative 
excess mass m e ). The tests are performed for a periodic system of 16 3 hydrodynamic cells of size 
Ax = 1, with fluid density p — 1 and shear viscosity 77 = 1 (in arbitrary units), with the particle 
started at position x = r = 10, and follow the motion of the particles to a time T — 80. In 
the left panel of Fig. [2] we show the error (51) for four different values of the excess mass for 
a spring constant k = 0.01, which implies a small Reynolds number Re ss r = u mBX Ax/v ~ 0.02, 
where u max is the maximal speed of the particle. As expected, the figure shows clear second-order 
convergence for neutrally-buoyant particles (m e — 0) , and a transition from essentially second-order 
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(at larger At) to first-order (at smaller At) accuracy in the other cases. Notably, for the larger 
time steps, which are of more practical interest, we see second-order convergence. As expected, 
the transition from second to first order of accuracy occurs at a larger At for particles that are far 
from neutrally-buoyant, and for m e /m f = 10 we see first-order accuracy over a broader range of 
time step sizes. 



In order to compare the three choices for the the slip correction, (43), (49) and Su n+ i — 0, in 
the right panel of Fig. [2] we compare the error for the three choices for the case of a large particle 
excess mass m e = I0m f and a larger Reynolds number, Re « 2 for k = 1, and Re « 20 for k — 500. 
We see qualitatively similar behavior as for the small Reynolds number case k = 0.01. Compared 



to the simple choice 5u" + 2 = 0, we see a modest improvement in the error when we use (43) for 



the largest Reynolds number, and we see a small improvement when we use (49) for intermediate 
time step sizes. Note that for our choice of parameters the viscous CFL number is /3 = At and the 
advective CFL is a ps /3 Re. For large Reynolds numbers the time step is limited by the requirement 
a < 1. 



Since the majority of the tests presented here are at small Reynolds numbers, we use (49) instead 



of the more expensive (43). For several of the tests we have also tried Su n+ 2 = and observed 
similar results. 



B. Compressible versus Incompressible 

In this section we apply our scheme to a standard test for the coupling of spherical particle 
of hydrodynamic radius R H to a compressible (38] 04, 66-68] or incompressible [3U [SU [69] fluid 
solver. The velocity autocorrelation function (VACF) 

C(t) = {v.(0)v.(t)) = ±{v(0)-v(t)) : (52) 
of a single free Brownian particle diffusing through a periodic fluid is a non-trivial quantity that 
contains crucial information at both short and long times. The integral of the VACF determines 
the diffusion coefficient and gets contributions from three distinct stages. Firstly, at molecular 
times equipartition dictates that C(0) = k B T/m, an important signature of fluctuation-dissipation 
balance that has proven challenging for several fluid-particle coupling methods |44[ [54"[ |66[ [69] . We 
recall that for our particle the effective particle mass m = m e + m f includes the mass of the fluid 
dragged with the particle m f , as well as the excess mass m e . The compressible inertial coupling 
method is able to reproduce the intercept k B T/m very accurately even for relatively large sound 
CFL numbers 138 . 
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Figure 3: Velocity autocorrelation function (VACF) (52) of a single particle with excess mass m e — mf 



normalized by ksT/m so that it should be unity at the origin for a compressible fluid. (Left panel) Com- 
parison between a compressible fluid for several different speeds of sound c (compressibilities), as well as an 
incompressible fluid (c — > oo). Vertical lines indicate the sound time scale t c = 2Rh/c, and the asymptotic 
power-law tail (t/ty)^/ 2 is emphasized in the inset, where t w = R 2 H /v is the viscous time scale. The tail 
matches the theoretical predictions for a rigid sphere with the same effective mass immersed in an incom- 
pressible fluid [S3]- All runs use a small time step size so the dynamics is well-resolved at short times. (Right 
panel) Comparison between different Schmidt numbers S c = v/x f° r a- n incompressible fluid. A deterministic 
calculation, corresponding to the limit S c — > oo, is also shown. In the legend, the time step size is expressed 
in terms of the advective CFL number a. The scaling of the time axes is adjusted to overlap the power-law 
tails (see text). 



On the time scale of sound waves, t < t c = 2R H /c, the major effect of compressibility is that 
sound waves carry away a fraction of the particle momentum with the sound speed c. The VACF 
quickly decays from its initial value to C(t c ) w k B T/m, where to = m e + dm s / (d — I) includes an 
"added mass" m f /(d— 1) that comes from the fluid around the particle that has to move with the 
particle because of incompressibility [501 E3 EO] • The initial decay of the VACF due to sound waves 
will appear to be instantaneous (discontinuous) if one increases the speed of sound to infinity. The 
incompressible inertial coupling method should produce an intercept C(0 + ) = k B T/m in agreement 
with our derivation in Appendix [A] [see Eq. (21)], rather than the equipartition result valid for 
a compressible fluid. For example, for m e = m s and d = 3, we expect C(0 + ) = 0.8k B T/m. This 
illustrates the subtle complexity of coupling a fluctuating fluid solver to immersed particles, even 
in the absence of external and interparticle forces. In the left panel of Fig. [3] we show numerical 
results for the VACF for several different speeds of sound, obtained using the algorithm summarized 
in Appendix ID! The approach to the incompressible limit c — ¥ oo is evident in the figure, and it is 
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Parameter 


Fixed S c runs 


Variable S c runs 


excess mass m e 


m e = rrif 


m e = mf 


grid spacing Ax 


10.543 


1 


grid size N 


41 3 


32 3 


fluid density p 


1 


1 


shear viscosity r\ 


0.5 


variable 


bulk viscosity ( 


0.5 


not relevant 


speed of sound c 


1 — oo 


oo 


Temperature fc^T 


1 


0.1 


viscous CFL ft 


10~ 5 - 10~ 3 


P = a^S c / (6tt) 


Sound CFL a s 


0.05-0.1 


not relevant 


Schmidt number S c 


48.2 


variable 



Table I: Parameters used in the compressible and incompressible simulations shown in the left panel of Fig. 
[3] (middle column), as well as the incompressible simulations shown in the right panel of Fig. [3] and in Fig. 
[1] (right column) . 

clear that our incompressible inertia coupling method correctly reproduces the limiting behavior 
(without however suffering from the severe time step limitation of compressible flow solvers). The 
parameters for these simulations are shown in Table [TJ 

At the viscous time scale, t > t v = pR 2 H /r], conservation of momentum (hydrodynamics) in the 
fluid introduces a memory in the motion of the particle and the VACF decays with a well-known 
asymptotic power-law tail ~ (t/t v )~ d/2 [71]. Any numerical method that solves the time-dependent 
Navier-Stokes equation (as opposed to the steady or time-independent Stokes equation, as Brownian 
or Stokesian dynamics do) ought to reproduce this power-law decay. The amplitude of the decay 
depends on the shape of the particle and is well-known for the case of a rigid sphere with stick 
boundaries [53]. We expect that the VACF for our blob particle will have the same value at the 
origin and the same power-law tail as a rigid sphere of radius R H and the same ratio m e /m f . 
For the "equivalent" rigid sphere we take m f = pV s where V s is the volume of the sphere, and 
m e = (p s — p) V s , where p 3 is the density of the sphere. The exact shape of C(t) will in general be 
different between a blob and a rigid sphere. In the inset of the left panel of Fig. [3] we compare 
the long-time behavior of the VACF of a blob particle to that of a rigid sphere, and see the same 
power-law behavior. Note that the rigid sphere theory shown here does not account for finite-size 
effects. In these tests we use a small time step size in order to study the properties of the spatial 
discretization in the absence of temporal truncation errors. 
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C. Diffusive Behavior 

We now investigate in more detail the long-time behavior of the particle velocity and displace- 
ment for larger time step sizes. Large At here means that at least one of the advective (a) or the 



viscous (/?) CFL numbers defined in (50) becomes comparable to unity. For a system at thermo- 



dynamic equilibrium at rest the typical value of the advection velocity to be used in (50) is the 
equilibrium magnitude of the thermal velocity fluctuations, 

y P Ax 3 

At long times, the motion of the particle is diffusive with a diffusion coefficient predicted by the 
Stokes-Einstein relation to be 

X * XSE = (tk (53) 
in three dimensions, where we recall that for the particular spatial discretization we employ R h ~ Ax 
in three dimensions. The exact coefficient depends on the system size, and for the system size we 
use here an excellent approximation is R h « Ax. 



We define here a dimensionless Schmidt number S c based on the prediction ( 53 ) , 



v _ v _ GnrfRn _ 67r/? 2 

X C Xse pk B T a 2 
rather than the actual diffusion coefficient x which in our model is a product of rather than an 

input to the calculation. The Schmidt number is an important quantity that measures how fast 

momentum diffuses relative to the particles. In many cases of interest S c 3> 1, which means that 



the dynamics of the particles approaches the Brownian (overdamped) limit (22), as discussed in 

Section |ITE~2"1 Note that the limit S c — > oo is the same as the deterministic limit k B T — > 0, in which 

fluctuations become a very weak perturbation to the deterministic dynamics. Another important 

dimensionless number is the "thermal" Peclet number 

VR f 



Xse 

which measures the relative importance of advection by the thermal velocity fluctuations to dif- 
fusion. We see that Pe ~ S\ /2 is directly related to the Schmidt number. A similar calculation 
also shows that the "thermal" Reynolds number is Re ~ S~ 1/2 . Therefore S c is the only relevant 
dimensionless number for a particle diffusing in a fluid at thermodynamic equilibrium. 

It is important to test how well our algorithm works for a range of Schmidt numbers. We 
expect the case of small S c to be the most difficult in terms of accuracy, since the particle can 
move a substantial distance (compared to the grid spacing) during a single time step, a — O(l), 
and the thermal and cell Reynolds number is also r = O(l). The case of large S c , on the other 
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hand, is the most demanding in terms of computational effort since particles barely move during 
a single time step, a <C 1, and 0(S C ) fluid time steps may be required to reach the diffusive time 
scale for (3 — 0(1). In order to investigate the long-time behavior we try to maximize the time 
step, but always keeping a < 1, specifically, here we set a = 0.25. For the largest Schmidt number 
we investigate, S c ~ 200, this value of a corresponds to a relatively large (3 w 0.81, and therefore 
we also try a smaller time step, corresponding to a = 0.1 and (3 ~ 0.33. The temporal integrator 
developed here cannot be used for /3 > 1 because the Crank-Nicolson temporal integrator we use 
for the velocity equation does not accurately resolve the dynamics of the small wavelength fluid 
modes |2SIH2]. 

1. VACF 

In the right panel of Fig. [3] we show the VACF for several viscosities and thus Schmidt numbers. 
The parameters for the runs are given in Table [TJ The standard theory for the tail of the VACF 
(long-time behavior) [53] implicitly assumes that S c ~> 1, and leads to the conclusion that for an 
isolated particle in an infinite fluid asymptotically C{t) » {t/t v )~ d/2 ~ {vt)~ m . A more complete 
self-consistent mode coupling theory [72] corrects this to account for the fact that while momentum 
diffuses around the particle the particle itself diffuses, and predicts that C(t) ~ [{x + v)t]~ d/2 [71] . 
This means that we expect the tails of the VACFs for different S c values to collapse on one master 
curve if we plot them as a function not of (i/t„) but rather of (1 + S^ 1 ) (t/t v ). This is confirmed in 
the right panel of Fig. [3] 

It is evident in Fig. [3] that these fluctuating calculations lead to noisy results for the tail of 
the VACF, making it difficult to see the behavior of the long-time behavior. Many researchers 
have chosen to calculate the VACF by performing a deterministic calculation, in which the particle 
is given a small initial kick in velocity, and then the deterministic algorithm is used to track the 
subsequent decay of the velocity. This is sometimes done because thermal fluctuations are not 
consistently included in the algorithm, or because the deterministic calculation is much faster and 
more accurate, not requiring as much statistical averaging. For linear systems, it is not hard to 
show that the correlation functions of the fluctuating system will be the same as the decaying 
correlations in the deterministic system. However, it is important to emphasize that this is not the 
case for nonlinear systems, and therefore, there is no a priori reason to believe it is strictly true 
for our nonlinearly coupled particle-fluid equations. Notably, the no-slip condition q = J (q) v is 
nonlinear because of the dependence of the local averaging operation on the position of the particle. 
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In the right panel of Fig. [3] we show the VACF obtained from a deterministic test, which can 
be thought of as the VACF in the limit of vanishing fluctuations, k B T — > (equivalently, S c —> oo). 
Note that in the deterministic test the magnitude of the initial velocity u Q of the particle has to be 
chosen to match the thermal kick, ||tt || 2 = dk B T/rh, which in practice means that the deterministic 
VACF has to be scaled so that it agrees with statistical mechanics at the origin. Due to the slight 
anisotropy and imperfect translational invariance of the spatial discretization, in principle even the 
deterministic result should be averaged over many initial positions and velocities of the particle. 
The VACF in the limit S c — > oo shown in Fig. [3] matches the fluctuating runs for the larger Schmidt 
numbers. Due to the lack of noise, it also clearly shows the long-time exponential decay in the 
VACF at times t ~ L 2 jv |54j . where finite-size effects become important. 



2. MSD 



In the left panel of Fig. [4] we show the time-dependent diffusion coefficient x(t) = Aq 2 (t)/ (2t) 
for a single particle in a periodic domain, where Aq 2 (t) is the mean square displacement (MSD) 

Aq 2 (t) = ([q x (t) ~ qM} 2 ) = ~ (\\Q(t) <?(0)|| 2 ) . (54) 
At long times, t >• L 2 /v, where L is the length of the periodic box, the VACF is expected to switch 
from power-law to exponential decay |54j . and the MSD should approach a diffusive behavior 
Aq 2 (t) R3 2\t. The long-time diffusion coefficient \ = limt->oo x(*) = f^ C(t)dt can be estimated 
most effectively by numerically integrating the discrete VACF up to some time t L , 

C(t)dt (55) 
which is also equivalent to numerically estimating the slope of Aq 2 (t)/2 at time t L . The value of 
t L should be chosen so as to balance statistical errors with the bias introduced by truncating the 
integral. Based on a statistical error analysis we have chosen to use t L = L 2 / (4z/) « 250 t v for our 
simulations. Note that for a finite system the integral of C(t) is well-defined in both two and three 
dimensions, even though the value of \ has a finite-size correction that diverges as bxL in two 
dimensions, and converges as in three dimensions |73j . 

For sufficiently large S c in three dimensions we expect \ to be well- approximated by the Stokes- 



Einstein prediction (53), independently of the inertia of the particle [53J. Obtaining xse requires 
estimating the effective hydrodynamic radius R H . We do this by turning off fluctuations, and then 
measuring the deterministic velocity of the particle under a small applied force in a periodic box, 
with an equal but opposite force applied homogeneously to the fluid to keep the center-of-mass 
from moving. Note that there are small variations in the drag as the particle moves through the 
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Figure 4: (Left panel) The ratio of the time-dependent diffusion coefficient x(i) = Aq 2 (t)/ (2t) to the 



long-time diffusion coefficient estimate xse given by the Stokes- Einstein formula ( 53 ) . No attempt is made 
to collapse all of the curves. The vertical line shows the time tx, = L 2 / (Av) at which we estimate the long- 
time diffusion coefficient x ~ X (Right panel) Comparison between the numerical estimate (symbols) 
for the long-time diffusion coefficient and the Stokes-Einstein prediction for several Schmidt numbers. The 



predictions of the self-consistent theory ( 58 ) are shown for comparison. The predictions of linear response 



theory (56) based on direct measurements of the mobility are shown for several Schmidt numbers. 



grid because the mobility — JC X S is not strictly translationally invariant, as discussed in more 



detail in Section IIIA1 and also Ref. [3]. Here we use an average value of the mobility over all 



positions of the particle relative to the grid. Note that for a finite system R H depends on the 
system size L due to the hydrodynamic interactions of the particle with its periodic images. This 
dependence is very- well captured by a simple theory [3J EH] , 



,R1 



R1 



Rh (L) = \l - 2.84-^ , „„. 
where R^ = (0.91 ± 0.01) Aa; is the numerical infinite system value 

In the right panel of Fig. [4] we show numerical estimates of the long-time diffusion coefficient as 



a function of the Schmidt number, obtained using a discrete form of (55). We observe deviations 



from (53) for moderate and small Schmidt numbers, and we observe good agreement for large S c . 
The diffusion coefficient obtained from the integral of the deterministic VACF (the limit S c — > oo) 



is within a couple of percent of the Stokes-Einstein prediction. The deviations from ( 53 ) for small 



S c are, at first glance, surprising. However, it should be noted that, with finite-size corrections 



taken into account, the Stokes-Einstein relation, as written in (53), can only be rigorously justified 



in the overdamped limit S c —> oo. More precisely, the separation of time scales required for the 
overdamped limit to hold rigorously is S c » (L/R H ) 2 > 1. 
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On the other hand, the Einstein formula, which is a linear response relationship between the 
diffusion coefficient and the mobility, x = {^bT) fi, is quite general, and follows from straightforward 
statistical mechanics arguments [73]. However, one must be careful here in defining the mobility /x. 
Mobility should be defined here as the coefficient of proportionality between a small applied force 
F on the particle (with the opposite force applied uniformly to the fluid to prevent center-of-mass 
motion) and the mean velocity of the particle, 

M = JimM (56) 
where the average here is taken over time (or, equivalently, over a nonequilibrium but steady-state 
ensemble that is a weak perturbation of the equilibrium ensemble). 

Because of the nonlinearity of the fluid-particle equations, thermal fluctuations can affect av- 
erage values. Therefore, we should not expect, in general, that the mobility measured at finite 
temperature is the same as the deterministic (k B T = 0) mobility. Indeed, in the right panel of Fig. 
[4] we show the predictions of (56) with (u) measured numerically under a small applied force (to 
ensure the system remains in the linear response regime). These predictions are in agreement with 
the direct measurements of the diffusion coefficient from the integral of the velocity autocorrelation 
function, as predicted by linear response theory. We therefore observe no violation of any known 
physical principles, so long as we recognize the fact that fluctuations cannot just be turned off. 
More precisely, at finite value of k B T the fluctuations do not necessarily average out, and can, in 
fact, affect the dynamics of macroscopic or deterministic observables. 

We are not aware of any theory that successfully predicts the form of the corrections to the simple 



Stokes-Einstein prediction ( 53 ) for moderate values of S c (see Ref. [75J for one attempt). The scaling 
of the tail of the VACF suggests a self-consistent equation of the form x{ v + X) = k B T j (6npR H ). 
In Ref. {73] a more systematic self-consistent theory is presented that suggests that the diffusion 



coefficient of the particle can be approximated as 



which is the solution to the self-consistent equation 



k P T ^ 1/2 



3npi>' 2 R H 



(1 + 25-^-1 , (57) 



In the right panel of Fig. [4] we observe agreement of numerical data with the self-consistent theory 
to within the statistical error bars. 

While the above results elucidate some important fundamental questions about the physical im- 
portance of fluctuations, it should be observed that small S c implies that the size of the blob particle 
becomes comparable to the length scale at which hydrodynamics breaks down (molecular diameter 
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grid spacing Aa; 


1 


grid size 


32 3 


fluid density p 


1 


shear viscosity 77 


1 


time step size At 


1 


temperature kgT 


0.001 


LJ strength e 


0.001 


LJ cutoff radius a 


2 


number of particles N 


1000 



grid spacing Aa; 


1 


grid size 


128 3 


fluid density 


1 


viscosity 


1 


advective CFL a 


0.01, 0.1, or 0.25 


viscous CFL f3 


9.2, 92, or 230 


excess mass m e 


m f 



Table II: (Left) Parameters used in the RDF simulations shown in the left panel of Fig. [5] (Right) Simulation 
parameters for the hydrodynamic interaction simulations presented in the right panel of Fig. [5j 

for liquids, or mean free path for gases). The physical fidelity of the fluctuating hydrodynamics 
equations used here should therefore be questioned at small S c . Nevertheless, our model allows 
us to separate hydrodynamic from molecular effects, and to test the predictions of self-consistent 
renormalization theory in two dimensions, where fluctuations play an even more prominent role 
[73]. 

D. Radial Distribution Function 



One of the most important requirements on any scheme that couples fluctuating hydrodynamics 
to immersed particles is to reproduce the Gibbs-Boltzmann distribution at thermodynamic equilib- 
rium. In particular, the probability distribution of the positions Q = {q 17 q 2 , . . . , q N } of a collection 
of N particles interacting with a conservative potential U (Q) should be 

P(Q)~exp[-^|>], (59) 
independent of any dynamical parameters such as viscosity or particle inertia. This follows from 
the balance between the dissipative and stochastic forcing terms and requires consistently including 
thermal fluctuations in the momentum equation. 

We verify that our incompressible inertial coupling algorithm gives the correct equilibrium 
distribution P (Q) by computing the radial (pair) distribution function (RDF) g(r) for a collection 
of colloidal particles interacting with a pairwise potential V(r), U (Q) = ^2a=i V ~ We 
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Figure 5: (Left panel) The equilibrium radial distribution function g(r) for a suspension of particles 



interacting with WCA potential (601 with cutoff radius a. The results for two different particle inertias are 



compared to a Monte Carlo sampling of the Gibbs-Boltzmann distribution. (Right panel) Hydrodynamic 
interaction force between two particles as a function of the interparticle distance. For large distances the 
Stokes mobility is recovered, and for moderate distances the next-order correction (Rotne-Prager mobility) 
is recovered, independent of dynamical parameters. At close distances a sharp (but not divergent) increase 
in the hydrodynamic force is observed. 



use the purely repulsive truncated Lennard-Jones (WCA) potential 



4c 



r < 2 1 / 6 a 



V(r) = { ^ ' (60) 
0, r > 2 1 < /6 cr 

In Fig. [5] we compare g(r) between a simulation where the particles are immersed in an incom- 
pressible viscous solvent, and a standard computation of the equilibrium RDF using a Monte Carlo 



algorithm to sample the equilibrium distribution (59). The parameters for these simulations are 



given in Table [TTJ For both m e = and m e = m f we obtain excellent agreement with the Monte 
Carlo calculations, even for the rather large time step size ft = 1. 



E. Hydrodynamic Interactions 

In this section we investigate the hydrodynamic interaction force between two particles in the 
deterministic setting, as done for a compressible fluid in Ref. [38]. In this test, we apply a force 
F on one particle toward the other particle, and the opposite force — F on the other particle, so 
that the center of mass remains at rest. The applied force is weak so that the Reynolds number 
Re < 10~ 3 and the flow is in the Stokes (steady-state) limit. As the particles approach, we measure 
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the relative speed of approach v r and compare it to the prediction of Stokes's law, 

= . (61) 

In the right panel of Fig. [5] we show results for F/F Sto]iCS as a function of the interparticle 
distance I. The simulation parameters are reported in Table [TTJ We performed several simulation 
to verify that the results presented here are independent of the excess mass and time step size. 
In the figure we compare the results of our calculations to a theoretical calculation based on the 
Rotne-Prager tensor for the mobility of a pair of particles in a periodic system |76j . We see that 
the Rotne-Prager correction correctly captures the behavior for the blobs for distances d > 3. 
At short distances there is a strong repulsive force between the blobs similar to the well-known 
"lubrication force" that develops as an incompressible fluid is squeezed out between two approaching 
rigid spheres. However, unlike the lubrication force between two rigid spheres, the hydrodynamic 
interaction force between two blobs does not diverge like (d — 2R H )~ 1 . This is expected because 
blobs do not have a well-defined surface; however, they are not point particles either, and they do 
squeeze the incompressible fluid in-between them as they approach each other. 



F. Large Reynolds Number 

As discussed in the Introduction, we believe that the blob particle model can be successfully 
used for simulations of particle-laden flows at moderate and large particle Reynolds number. In the 
Re P — > limit, we showed that the perturbative flow created by the blob particle agrees with the 
Stokes solution for the steady flow around a no-slip sphere of hydrodynamic radius R H = 0.91 h for 
the three point kernel and a staggered grid solver [38]. This is consistent with other calculations 
utilizing a four-point kernel and a non-staggered grid |61j . However it is not a priori clear how 
consistently the blob hydrodynamic behavior will be at large Re P . The drag force provides a non- 
trivial test because it captures the average effect of the perturbative flow. In a previous study [38] 
we observed that the drag on the blob is consistent with that of a rigid sphere up to Re P < 10. 
This study was performed using a compressible flow solver, and ensuring a low Mach number 
was prohibitively expensive at larger Re P . Using the incompressible (zero Mach number) solver 
developed here we now study larger Re P . 

In the top left panel of Fig. [6] we show the drag force on a blob particle in a periodic domain as 
a function of Re, normalized by the Stokes limit drag (Re — > 0). We estimate the drag coefficient by 
dragging the blob with a constant applied force and measuring the average velocity v = (it) along 
the direction of the applied force at long times. The particle interacts with its periodic images 
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Figure 6: ( Top left) Drag force on a blob particle in a periodic domain as a function of the particle Reynolds 
number Rep = 2Rh (u) /is, normalized by the Stokes drag (Rep — > 0) (symbols). Results for the incom- 
pressible and compressible solvers (see Appendix [D] and Ref . [38] ) are compared with the empirical law for 
rigid no-slip spheres (line) [77] . The size of the domain box in cells is indicated in the legend. (Top right) 
Out-of-plane vorticity iso-contours at Rep = 1.5 for a box of 32 3 cells, and (Bottom left) at Rep = 137 in 
a long box of 512 x 32 x 32 cells. (Bottom right) A snapshot of an unsteady (nearly oscillatory) flow in a 
box of 32 3 cells at Rep = 70 where the particle sheds vortices due to the interaction with the wake from its 
image. The simulation parameters are as in Table \U\ but m e = and advective CFL a = 0.2. Error bars 
are estimated to be less than 5%. 



approximately after time t l = L/v (where L is the box length in flow direction), while viscous 
transport around the blob requires a time longer than t v = R 2 H jv to settle down. In order to mimic 
the behavior of an isolated blob in an infinite medium, we must have t l /t v > 1, or, equivalently, 
L/R H > Re p . We therefore performed the calculations in a box of 2™ x 32 x 32 cells, with 2" > 3Re F ; 
the size of the simulation is indicated in the legend in the figure. In the top left panel of Fig. [6] we 
compare the numerical results with the empirical law for the drag on a rigid sphere with a no-slip 
surface [77]. The agreement is remarkably good over the studied range < Re P < 324. 
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Encouraged by this result we now briefly analyze the structure of the perturbative flow around 
the blob, and defer a more detailed study for future work. In the remaining panels in Fig. [6] we 
show iso-contours of the vorticity perpendicular to the snapshots' plane at a few values of Re P . For 
Re P ~ the fully symmetric pattern observed around the blob agrees with the Stokes solution at 
distances larger than 2R H [38]. The fore-and-aft symmetry of the Stokes flow becomes distorted 
by advection at moderate Reynolds Re P ~ 1 leading to the so-called Oseen flow |16j , as observed in 
the top right panel of FigJBJfor Re P = 1.5. For 20 < Re P < 270 a transition with symmetry breaking 
leads to a steady axisymmetric "double-thread" structure [751 E2] • Although we have not studied 
the transition in detail, our blob model qualitatively reproduces a steady axisymmetric wake with 
a bifid vortex trail, as illustrated in the bottom left panel of Fig. [6] for Rc P = 137. This type of 
wake is observed above Re P > 10 and its topology is maintained at least up to Re P ~ 300. For flow 
around a rigid sphere a Hopf transition to oscillatory flow leading to vortex shedding takes place at 
Re P « 270 \78\ 179]. For the blob we observe small oscillations of the wake, without vortex shedding 
at Re P > 300, indicating that some possible transition to unsteady flow could be induced by a small 
perturbation. In fact, in particle-laden flows typically Re P < 100 [32j, and the induction of vortex 
shedding due to perturbations from the wake of other particle is a more relevant vorticity source. 
Using smaller boxes, where the particle interacts with its image, we frequently observed induced 
vortex shedding for Re P > 70, as illustrated in the bottom right panel of Fig. [6j 

It is a remarkable fact that the ICM blob minimal-resolution model can produce wakes con- 
taining many of the features of realistic flows. The thickness of the viscous (Oseen) layer around 
a sphere decreases like i?/Re P [IB] : therefore, this layer is completely unresolved for Re P > 1 in our 
model. However, the "local" no-slip constraint is able to capture the non-linear velocity-pressure 
coupling which dominates the drag at large Re P . A somewhat similar scenario was observed in sim- 
ulations of ultrasound-particle interaction [38] . In the inviscid regime (sound frequency faster than 
is/R 2 H ), an excellent agreement with the theory was observed even in cases where the sound- viscous 
boundary layer around the blob was unresolved. 

V. CONCLUSIONS 

In this work we described a bidirectional coupling between a point-like "blob" particle and an in- 
compressible fluctuating fluid, building on prior work by some of us in the compressible setting |38j. 
At the continuum level, the proposed model includes inertial and stochastic effects in a consistent 
manner, ensuring fluctuation-dissipation balance and independence of equilibrium thermodynamic 
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properties on dynamical parameters. We constructed a second-order spatio-temporal discretization 
that tries to preserve the properties of the continuum model as well as possible. 

Through numerical experiments, we showed that the proposed inertial coupling method (ICM) 
can consistently describe the dynamics of "blob" particles in fluid flow over a broad range of particle 
and fluid Reynolds numbers, from Brownian motion to convective regimes. We demonstrated 
that the method reproduces well-known non-trivial effects of particle inertia on the short-time 
dynamics of Brownian particles, while also capturing the long-time (Brownian) diffusive dynamics 
and the associated equilibrium distribution correctly. Remarkably, we found the minimally-resolved 
blob model to reproduce non-trivial features of turbulent flow around a rigid sphere, including 
the non-Stokes drag at moderate Reynolds numbers and the transition to turbulence at higher 
Reynolds numbers. As such we believe that the method presented here can be applied to model 
the dynamics of dilute and semi-dilute colloidal suspensions and polymeric fluids over a broader 
range of conditions than existing methods [3l [20 1211 ESI EZ] • 

The algorithm we described here was specifically optimized for periodic boundary conditions. In 
particular, we constructed a semi-implicit temporal integrator that only requires a few applications 
of the FFT algorithm per time step. This enabled us to implement the algorithm on GPU platforms, 
achieving excellent performance with little development effort. However, this simplicity was not 
without a cost. Firstly, we had to make several approximations in order to avoid iterative linear 
solvers and use a fixed number of FFTs per time step. Notably, we had to assume that particles 
were far away from other particles compared to their size (dilute conditions). Secondly, in many 
problems of interest the fluid is confined in non-periodic geometries such as channels and periodic 
boundary conditions are not appropriate. It is not difficult, at least in principle, to adopt our 
algorithm to non-periodic geometries and to dense collections of blobs. This requires developing 



specialized preconditioned Krylov linear solvers for solving the "inertial" Stokes problem (34) in non- 
periodic geometries, similar to the fluid-only solver developed by one of us [371 [5B] and implemented 
in the IB AMR framework |57j . In future work we will explore whether the "splitting" technique 



developed in this paper can be used to construct an effective preconditioner for solving ( 34 ) using 
a Krylov method. 

The temporal integrator algorithm we developed in this work can accurately resolve the inertial, 
viscous, and fluctuating dynamics of particles immersed in an incompressible fluid if the viscous 
and advective CFL numbers are less than unity. Some further work is required to tackle the 
large separation of time scales present in many situations of practical interest. For example, in the 
Brownian dynamics limit there is an increasing separation of scales between the particle movement, 
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the viscous damping, and the inertial dynamics. For the second-order temporal integrator that we 
presented here, the case of large Schmidt number S c ^> 1 is the most demanding in terms of 
computational effort. This is because particles barely move during a single time step, a<l, and 
0(S C ) fluid time steps may be required to reach the diffusive time scale for /3 = 0(1). For the case of 
a neutrally-buoyant particle (m e = 0) and periodic boundary conditions Atzberger et al. I28\j have 
developed a specialized exponential integrator that tackles the large separation of time scales that 
arise when S c 3> 1. Future work will consider extending such techniques to non-periodic systems 
and to the case m e ^ 0. 

In ICM the fluid-particle coupling is expressed as a no-slip constraint equating the translational 
particle velocity with the local fluid velocity. This implies that only the monopole term (stokeslet) 
is included in the fluid-particle force. As a consequence the present approach can only accurately 
resolve the fluid flow at distances larger than the typical size of the particles. However, it is 
a remarkable fact that with such minimal resolution permits to capture so many hydrodynamic 



effects in a qualitatively, and, in some cases, quantitatively correct way (see Sec. IV and also 
Ref. [38]). Even though lubrication flows (at low Re P ) or viscous boundary layers (at large 
Re P and/or high forcing rates) are unresolved, the locally averaged no-slip constraint proves to 
be remarkably robust in capturing their essential hydrodynamics and permits to go beyond the 
Stokes limit. Crucial to this success is the fact that we employ a carefully-constructed spatial 
discretization, combining nearly grid-invariant Lagrangian-Eulerian transformations, local energy 
and momentum conservation, and a staggered discretization of the incompressibility condition and 
the fluctuating stresses. 

It is not difficult to extend our approach to also include the anti-symmetric component of the 
dipole (rotlet) stress [29J. Firstly, particle rotational degrees of freedom would need to be added to 
the blob description, along with an angular velocity a; and an associated excess moment of inertia. 
We would need to impose an additional rotational no-slip constraint, requiring that the particle 
rotate with the locally-averaged angular velocity of the fluid, 

u,= l -J{V xv), (62) 
and distribute the fluid-particle torque t (Lagrange multiplier) as a force density f T = - V x (St) /2 
in fluid equation. This type of approach has already been employed in the deterministic context to 
model suspensions of neutrally-buoyant semi-rigid rods |8U| . To our knowledge, such an approach 
has not yet been applied to fluctuating hydrodynamics and the resulting rotational diffusion has 
not been investigated. 
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An additional rigidity constraint would be required to also constrain the locally-averaged de- 
formation tensor, and thus consistently include the symmetric components of the dipole (stresslet) 
force terms Unlike the smooth Gaussian kernels used in the force-coupling method [29J, the 
existing Peskin kernels do not have well-behaved derivatives. Therefore, it appears necessary to 
generalize the IBM kernels to enable the local averaging ("interpolation") of spatial derivatives in a 
reasonably translationally-invariant manner. The inclusion of thermal fluctuations of the stresslet 
requires careful consideration even at the continuum level, and will be the subject of future research. 

In a different spirit, the approach used here for single "blobs" can be extended to account for 
the finite extent and shape of arbitrarily-shaped rigid bodies immersed in fluid flow. One approach 
that has been successfully employed in the deterministic setting is to construct the immersed body 
out of a collection of blobs constrained to move rigidly [81] . Future work will consider the inclusion 
of thermal fluctuations in this type of approach, and the minimal amount of resolution required to 
capture the geometry of immersed particles. 

At the other extreme, our method can be used to provide a coarse-grained model for hydrody- 
namics at very small scales, for example, for the Brownian motion of a small molecule suspended 
in a simple solvent. At such small scales, the suspended particles do not have a fixed, or even 
a well-defined shape, and there are many competing and sometimes canceling effects: normal 
and tangential slip [82], breakdown of Navier-Stokes hydrodynamics, non-Gaussian fluctuations, 
non-Markovian effects, etc. Our simple fluid-particle coupling model can be used to isolate hydro- 
dynamic from non-hydrodynamic effects and study basic physics questions about the importance 
of inertia and fluctuations on Brownian motion, going beyond the uncontrolled approximations 



required by existing theoretical approaches. Notably, our results in Section IV C showed that for 
small Schmidt numbers nonlinear effects become important and lead to a non-trivial contribution 
of the thermal fluctuations to the mean fluid-particle force. Specifically, the mobility of a particle 
in a fluctuating fluid was found to differ from that in a deterministic (Stokes) fluid, thus leading 
to a deviation of the diffusion coefficient from the standard Stokes-Einstein prediction. More care- 
ful investigations are required to assess how well self-consistent mode-coupling theories can model 
this effect, as well as to study the influence of particle inertia (density contrast), random slip (see 
Appendix IB]) , and spatial dimensionality. 
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Appendix A: Periodic Boundary Conditions 



The pressure- free fluid-only equation (17) can formally be written in a form suitable for direct 
application of standard numerical methods for integrating initial value problems, 

d t v = P J(Vf + VSF), (Al) 
although this form is only useful if one can actually compute the action of the inverse inertia tensor 
p~g. It turns out that this is possible for periodic systems. 

To see this, we expand p^ 1 into a formal series, 

pj = p- 1 \l - m eP - x VSJV + {m e p- x VSJVf - {m e p- x VSJVf + ...]. (A2) 
Observe that {VSJV) n = VS {JVS) n JV involves the n-th power of the dx d matrix AV _1 = JVS, 
where we made use of the fact that V 2 = V . Recall from ([!]) that JS = AV' 1 is a scalar constant 
related to the inverse of the particle volume, which is independent of the position of the particle 
for the particular kernel function used herein. In principle, the matrix AV could depend on the 
position of the particle because of the appearance of P, which implicitly encodes the boundary 
conditions. However, periodic systems are translationally invariant and therefore AV cannot 
depend on the position of the particle and is simply a constant d x d matrix. 

By performing a relatively-straightforward calculation in Fourier space it is possible to show 
that for periodic systems much larger than the kernel extent AV is simply a multiple of the 
identity matrix, 

AV _1 = JVS= ^— -AV' 1 !. (A3) 
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The prefactor (d—l)/d in the relation accounts for the elimination of the longitudinal (compressible) 
velocity mode by the projection operator. We can therefore define the effective particle volume 
accounting for incompressibility to be 



AV = - 

a - 1 

and accordingly, define the effective mass of the fluid dragged with the particle to be rh f = pAV = 
dm,f/ (d— 1), and the effective particle inertia to be m = fhf + m e . In three dimensions, the added 
fluid mass due to incompressibility is to//2, which is a well-known result for rigid spheres immersed 
in an incompressible fluid [S3 ED] • 

Note that for periodic systems where some of the dimensions of the unit cell are comparable to 



the kernel width (A3) is only an approximation. Notably, for a three-dimensional periodic box of 



shape L x x L y x L z , if L z <c a the value of AV converges to the two-dimensional result, 



AV~v = AV' 1 



\ o o 



i 



1 

rather than the three-dimensional AV 3D = (2AV _1 /3) I. 



(A4) 



Using ( A3 ) it is possible to simplify the infinite series ( A2 ) and obtain the important result for 

(A5) 



a single particle immersed in a periodic fluid, 

{ P i + 1 n e rsjry 1 



pj 



m e AV 



VSJVj 



Using this result we can rewrite (Al) in the simple form 

m e AV 



pd t v = V\I 



SJ Tf 



m 



(A6) 



which is useful for analysis and for numerical approximations. 

It is important to point out, however, that when there are many particles present, there is no 



simple formula for p cti . That is, we cannot just add a summation over all particles in (A5). This is 



because there are cross-terms between two particles i and j in the infinite series ( A2 ) involving the 
operator J{PSj. If the particles are not overlapping, meaning that the kernel functions of the differ- 
ent particles have disjoint support, then JiSj = 0, however, this is not true for incompressible flow 
because the projection V is a non-local operator involving the inverse Laplacian. Nevertheless, the 
cross terms decay fast as particles become well-separated from each other. Specifically, theoretical 
calculations suggest that to leading order J{PSj decays with the distance between the two particles 
like dipole and quadrupole terms, and is thus expected to be very small in semi-dilute suspensions 
|52j . In many problems of practical interest there are repulsive forces between the particles that 
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will keep them from coming close to each other, and the approximation J{PSj ~ (^^) ^ ^> 

' ^ 771 ■ 



PeS 



(A7) 



will be quite accurate even for multiparticle systems. We investigate the accuracy of the approxi- 



mation J iPSj = for i 7^ j numerically in Section III A 1 



Appendix B: Fluctuation-Dissipation Balance 



In this Appendix we demonstrate that the coupled fluid-particle equations written in the form 



(17) 



d t v 



dq 

~dl 



PlxT 
r/V 2 v + V 



p(v • V) + m e SJ v ■ — J 
{k B Tri)^ (W + VW T )]} 



v + SF 



= Jv, 



(Bl) 
(B2) 



obey fluctuation-dissipation balance with respect to the Gibbs-Boltzmann distribution with coarse- 
grained free energy 



H (v, q) = - j v T p cff v dr + U (q) , 



where the effective fluid inertia operator p cff is given in (18). The calculations here will follow the 



techniques described in detail in Ref. [51] (ignoring boundary terms), and are purely formal in 
the continuum (infinite-dimensional) setting. More precisely, the equations we write should really 
be interpreted as a short-hand notation for a spatially-discretized system in which there is a finite 
number of degrees of freedom [27] . 

It is well-known |47| l5"T] that within an isothermal and Markovian approximation the generic 
form of evolution equations for a set of macroscopic variables x has the Ito Langevin form 

d t x = -N (x) — + (2k B T) ^ B (aj) W(t) + (k B T) — ■ N* (x) , (B3) 
where W(t) denotes a collection of independent white noise processes, star denotes an adjoint, 
and (d x ■ N*) k = dN kj /dxj in indicial notation. We will suppress the explicit dependence on x and 
usually write the mobility operator as N = N (x). The fluctuation dissipation balance is contained 
in the relation 

BB* = M = -(N + N*), 



and the last term in (B3) is a "thermal" drift which ensures that the dynamics obeys 



detailed-balance (time-reversibility) at equilibrium with respect the Gibbs-Boltzmann distribution 

Z- 1 exp [-H (x) j k B T] . 
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In our case, the coarse-grained variables are x = (v,q) and the thermodynamic driving forces 
are given by the functional and partial derivatives 



dH 

- = -F{ q) + m e j[v.-J\ v . 



(B4) 
(B5) 



Note that the kinetic term 



,SJ (v- J^J)d appearing in (Bl) term follows from dH/dq, and 



therefore the approximation Oj w in ([6]) is not consistent in general. 

1. Mobility Operator 



The mobility operator that gives ( B1|B2 ) from the thermodynamic driving forces ( B4|B5 ) can 
be symbolically written in the block form 



N = 


pjN NS pJ 


pjrs 


, B = 


PeS B NS 




JVpJ 






B B o 



where the subscript NS denotes the corresponding operators for the fluctuating Navier-Stokes 
equations without any immersed particles [H]. In the form of the equations that we employ 
AT BD = B BU = 0. Fluctuation-dissipation balance then follows from the corresponding property 
of the pure fluid equations, B NS -B^ S = (7V NS + AT* S ) /2, since there is no dissipative terms in the 
particle equation. The fluid-particle coupling contribution is non-dissipative or skew-adjoint, which 
follows from the antisymmetry between to the lower left and upper right blocks in N. 

It is, however, consistent with the general framework of augmented Langevin descriptions and 
fluctuation-dissipation balance to allow for a nonzero B m with 7V BD = £? B d-B B d- This does not 
affect the fluid equation but changes the particle equation. For example, the simple choice 

B BD = ^/(I, N BD = CI, 
leads to a modified equation of motion for the particle, 



dq 

~dt 



v + CF(q) + y/2{k B TW q , 



where W q (t) is a white-noise "random slip". If v = the above equation is recognized as the usual 
equation of Brownian dynamics with friction coefficient £ -1 and mobility £. It is therefore expected 
that adding a random slip component would increase the diffusion coefficient of the particle by 
(k B T. This can be used to tune the diffusion coefficient to some target (e.g., experimental) value. 
In this work, we fix ( — 0, and in future work we will consider second-order algorithms for £ > 



54 



for the case of a neutrally buoyant particle, m e = 0. 



Thermal Drift 



Complications arise because of the presence of the last term in (B3) (proportional to k B T), 
which is non-zero because J and S depend on q. This "thermal drift" term arises because u is 
eliminated from the description as a variable adiabatically-slaved to v, neglecting the fact that the 
particle position fluctuates rapidly due to the fluctuations in u. This term was identified in Ref. 
|27j as missing from the formulation of the Stochastic Immersed Boundary Method [28], which 
is equivalent to our formulation with m e = 0. The missing contribution is an extra term in the 



velocity equation (Bl) of the form 



One can, in principle, numerically evaluate this term without requiring any derivatives by using 
"random finite-differences" via the identity 

~ ■ {JVp-J) = lim er 1 {[p-J (q + Aq) VS (q + Aq)] Aq - [pj (q) VS (q)} Aq) Aq . 
where the expectation value is with respect to the small random particle displacement Aq = eW q , 
and the components of W q are independent standard normal variates. 



For periodic boundaries, using (A5) we can simplify, for any q, 

rhf 
fn 

and therefore the thermal drift 



[Pelr 1 (?) VS (q)} Aq = ^-V [S (q)} Aq, 



/th = (k B T) lim e" 1 ([S (q + Aq) - S (q)] Aq) Aq = -^VV S (k B T) = 

vanishes as the projection of a gradient of a scalar. Since we only consider periodic systems, we 
have ignored the thermal drift term in this work. 

Note that some of the mathematical manipulations used above rely on continuum identities 
which are not necessarily obeyed by the discrete operators. In particular, for the spatial discretiza- 
tion that we employ, the discrete vector field w = -^S is not equal to a discrete gradient of a scalar. 
However, numerical observations suggest that it is very close to a discrete gradient of a scalar field, 
in the sense that application of the discrete projection reduces the magnitude by several orders, 
HPiull <C regardless of the position of the particle relative to the fluid grid. This is one more 
case in which we find that the Peskin discrete local averaging and local spreading operators closely 
mimic the properties of the continuum operators even though they were not specifically designed 
with the staggered-grid discrete projection V in mind. 
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Appendix C: Semi-Implicit Inertial Velocity Correction 



In this Appendix we explain how we solve the linear system ( 35|36|37|38 ) using only FFT-based 
linear subsolvers. 



First, we solve (35) for v n+1 and w n+ 2 using a projection algorithm and FFTs as explained 

(CI) 



above. Then, we eliminate u n+1 in (37) using the no-slip constraint (38) to obtain 

( 



At A" 



J ( v n+1 + A^+J ] + A«"+5 



Substituting this expression into (36), we get the fluid correction equation 

f AtV(A7r n+ 2) = m e S (u n - Jv n+1 



Av r ' 



Am" 



(C2) 



where L = V 2 denotes the discrete Laplacian operator. This is nothing but a rewriting of (31 ) and 
just as difficult to solve, so we need to make further approximations. 



We begin solving ( C2 ) by writing 



pi — 2~ ? ?-k + m cSJ = (pi + m e SJ) (i — 2~ 1/ I J 



m P vAt 



SJL. 



(C3) 



If we ignore the last term we obtain an approximation of order 0(m e At). We can improve the 
accuracy for viscous-dominated flows if we time lag the last term, that is, evaluate it during the 



previous time step. Using (C3) in (C2) along with time lagging we get the modified fluid correction 
equation 



(/>/ + ,„, SJ) [I--£uL } Ac" 3 -| A/ V ( Att ■'>-)= S 



mAu n - Jv n 



5u r ' 



= SAp, (C4) 



where we have time-lagged the last term in (C3) and denoted 



Note that if one approximates Su n+ i = (lowering the order of accuracy), then as At — > oo we 
get Av n+ i —> since both pi + m e SJ and I — AtvL/2 are positive semi-definite matrices. This is 
consistent with physical intuition that inertia should play a negligible role at long time scales, and 
implies that the algorithm will be unconditionally stable in the absence of advection. 

Denote Av = ( I - AtvL/2) Av n+ ^, and note that V ■ (Av) = V ■ (Av n+ i^j = with periodic BCs 



(C5) 



since the divergence and the Laplacian commute. Therefore, we can solve (C4) by first solving 

(pi + m e SJ) Av + AtV (A7r" + 3) = SAp, 
subject to V • (Av) — 0, and then solving 

fl - ^ vL\ Av n+ ? = Av, (C6) 

using FFTs. 



We solve (C5) using the pressure-elimination procedure described in Appendix [A] to obtain the 
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equivalent pressure-free equation 

(pi + m e VSJV) Av = VSAp, 
where V denotes the discrete projection operator [39]. By assuming the particles are sufficiently 



far away from each other, we can employ the approximation ( 26 ) to approximately solve ( C5 ) using 

-VS (u n - Jv n+1 -5u n+ ^ 



a the discrete equivalent of ( A6 ) , 



m e m f 



777 * 

A* = -^VSAp = 

pm p(rrif + m e ) 



(C7) 



where we recall that rh f = dm f / (d — 1). 

In order to ensure strict momentum conservation even in the presence of the approximation 



(26), we rewrite (C6) using (C5) in the form 



pi- ^*-r)L ) A» n+ i + A/ V (At" i ) = SlAp - JAv). 



subject to V • yAv" + 2 J = 0. This is consistent with (36) if we set 
A< (a™+2 +F n+ ^ =m 



J (¥' 



Av) +5u n+ 2 -u r ' 



(C8) 



(C9) 



which is to be compared to (CI). The linear system (C8) can be solved using the same techniques 



as used to solve (35). 



Having determined Ad" + 2 we can update the fluid velocity (momentum), v" +1 — v n+1 + Av" + i . 



If desired, we can update the particle momentum in a conservative manner by substituting (C9) 



into (37) to get 



t n+1 = J (v n+1 + Av) + Su n+ ^ . 



(CIO) 



From (CIO) we see that the approximations we used lead to a violation of the no-slip constraint 



(38) 



vAt 
~2~ 



which is of O (m e At 2 ) if Av is smooth in time, Av n+ ^ = A«""J + O(At). This means that our 
procedure approximates the solution of ( 35|36|37|38 ) with slip A«" +1 = O (m e At 2 ). Therefore, 
in the deterministic setting, for sufficiently smooth flows, the temporal integrator summarized in 



Section III C is expected to be second-order accurate. 



Appendix D: Compressible Algorithm 

In this Appendix we propose a modification of the first-order temporal integrator developed 
in Ref. [38] , following a similar approach to the incompressible algorithm summarized in Section 



IIIC The algorithm summarized below is still only first-order accurate if m e ^ 0, however, it is 
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second-order accurate for m e = unlike the algorithm used in Ref. [38J. Numerical results show 
that the modified algorithm below is a substantial improvement over that used previously [38j . 

1. Estimate the position of the particle at the midpoint, 

q n+l = q n + ^-J n V n , (Dl) 

and evaluate the external or interparticle forces F n+ ' . 

2. Solve the coupled density and unperturbed momentum equations [83J 

D tP = -p(V-t>) 
p(D t v) = -c 2 T Vp- V ■ <t + S n+h - F n+ ^- 
using a step of the third-order Runge-Kutta algorithm described in Ref. [39], to obtain the 
density p n+1 and the unperturbed fluid velocity v n+1 . During this step we treat the force 
density S" + 5f"+j in the momentum equation as a constant external forcing. 

3. If m e = 0, set v n+1 = v n+1 and skip to step[7| 

4. Calculate the momentum exchange between the particle and the fluid during the time step 



At 



+ = Ap = ™ emi m (j n+1 -v n+1 - u n ) , (D2) 



m e + m,f 

where the mass of the dragged fluid is estimated from the local density as 

m f = AV J n+ 2p™ +1 . 



5. Update the particle momentum, 

l l__^jn+i in +l_ u ny ( D3 ) 



u n+1 = u n + — = u n 



6. Update the fluid velocity to enforce the no-slip condition u n+1 — J n+2 v n+1 , 

v n+1 = v n+1 - — S n+ ?Ap = v n+1 + AVS n+ i (u n+1 - J n+ h n+1 ) (D4) 

rrif v J 

Note that this conserves the total momentum since 

pn+l^n+l _ p n+l _ _ (jn+^n+A" 1 f p n+l gn+ \ Ap dj , = _ Ap 



7. Update the particle position, 



? n + l _ q n + At jn+ 1 ^ n + 1 + ^ 



There are some differences between this algorithm and the incompressible algorithm summarized 



in Section III C Notably, viscosity is handled explicitly in the compressible algorithm to avoid 
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costly linear solvers. In the incompressible algorithm the no-slip condition is violated slightly 
while here the velocity correction v n+1 - v n+1 is designed to enforce the no-slip condition exactly 
even in the presence of density variations. Another important difference is that the compressible 
algorithm presented here is not second-order accurate even for small Reynolds number because 
viscous corrections of O (vAt) are not taken into account when computing X n+ ^. Since the time 
step is typically strongly limited by propagation of sound, typically |?« 1 and there is little need 
for higher-order handling of the viscous terms. 
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